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Abstract 

The  sharper  variation  points  of  a  signal,  measured  at  different  scales,  can  be  detected  from 
the  local  maxima  of  its  wavelet  transform.  We  describe  an  algorithm  that  reconstructs  one- 


dimensional  signals  and  images  from  their  sharper  variation  points  at  the  scales 


2' 


.  This 


algorithm  reconstructs  exactly  images  from  their  multiscale  edges.  It  provides  a  verification  of 
Marr's  conjecture  concerning  the  completeness  and  stability  of  multiscale  edges.  We  also  prove 
that  the  evolution  across  scales  of  the  wavelet  maxima  characterizes  the  local  shape  of  the  signal 
sharp  variations.  We  can  thus  not  only  detect  edges  but  also  classify  them.  The  wavelet  maxima 
representation  is  a  new  reorganization  of  the  image  information  that  enables  us  to  develop  algo- 
rithms uniquely  based  on  edges  for  solving  image  processing  problems.  We  discuss  more  particu- 
larly compact  image  coding,  texture  discrimination  and  pattern  recognition.  Although  we 
emphasize  image  processing,  this  representation  also  opens  new  approaches  to  the  processing  of 
many  types  of  one-dimensional  and  two-dimensional  transient  signals. 


This  research  was  supported  by  the  NSF  grant  IRI-890331  and  Airforce  grant  AFOSR -90-0040. 
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1.  Introduction 


Points  of  sharp  variations  are  often  among  the  most  important  features  for  analyzing  the 
properties  of  transient  signals  or  images.  In  images,  they  provide  the  locations  of  the  edges 
which  are  generally  used  for  describing  the  structural  information  of  object  boundaries.  Edge 
detection  is  generally  viewed  as  a  process  which  reduces  the  amount  of  data  representing  the  sig- 
nal but  which  does  not  preserve  the  whole  signal  information.  In  this  paper,  we  show  that  muliis- 
cale  edges  can  provide  a  complete  and  stable  description  of  signals.  The  model  is  based  on  ihc 
wavelet  formalization  of  multiscale  transforms.  This  characterization  has  many  applications  for 
pattern  analysis  but  also  provides  new  approaches  to  classical  signal  processing  problems  such  as 
compact  coding  and  enhancement.  We  concentrate  on  image  processing  applications  although 
the  model  is  described  for  both  one  and  two-dimensional  signals. 

The  measurement  of  a  signal  local  changes  must  be  defined  with  respect  to  a  reference 
scale.  The  scale  characterizes  the  neighborhood  size  where  theses  variations  are  computed.  In 
signals  which  include  meaningful  structures  of  very  different  sizes,  one  needs  to  modify  this  scale 
parameter.  At  large  scales,  we  measure  the  signal  variations  over  large  neighborhoods  whereas  at 
fine  scales  we  obtain  the  variations  due  to  the  firmer  structures  of  the  signal.  The  concept  of  mul- 
tiscale has  been  studied  within  the  computer  vision  community  but  also  in  mathematics,  physics 
and  signal  processing  in  general.  Recently,  this  effort  converged  in  a  more  coherent  mathemati- 
cal theory  called  the  wavelet  theory.  The  wavelet  model  formalizes  some  ideas  already  existing  in 
the  computer  vision  literature  [35]  but  also  provides  the  mathematical  tools  necessary  to  go 
beyond  the  ad-hoc  stage  which  often  limits  multiscale  algorithms. 

In  the  first  section,  we  summarize  briefly  the  most  relevant  matiiematical  results  of  Uie 
wavelet  model.  For  digital  computations,  the  wavelet  transform  must  be  discrctizcd.  The  com- 
pleteness and  inner  redundancy  of  linear  discrete  wavelet  transforms  are  now  well  understood 
[8,28].  However,  these  representations  are  not  translation  invariant.  When  a  pattern  is  translated 
in  die  image,  the  wavelet  coefficients  are  not  translated  but  totally  modified.  It  is  thus  difficult  to 
build  pattern  recognition  algoritiims  from  such  descriptors.  This  distortion  with  translation 
appear  in  all  the  multiscale  pyramidal  representations  used  in  computer  vision  such  as  the  Lapla- 
cian  pyramid  [2]  or  the  DOLP  transfonm  [4].  For  pattern  analysis,  we  cleariy  need  to  build 
descriptors  Uiat  translate  wiUi  the  signal.  We  introduce  an  adaptive  sampling  of  the  wavelet 
transform  based  on  its  local  maxima.  This  sampling  is  translation  invariant  and  for  a  particular 
class  of  wavelets,  the  maxima  indicate  Uie  position  of  sharper  variations  in  die  signal  amplitude. 
This  sampling  is  therefore  equivalent  to  a  multiscale  edge  detection.  In  this  model,  we  decom- 


pose the  signal  on  a  dyadic  sequence  of  scales 
tion. 
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g2  to  simplify  die  computer  implcmcnia- 
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An  important  issue  in  multiscale  analysis  is  to  relate  the  local  properties  of  the  signal  to  the 
evolution  of  the  transform  values  when  the  scale  varies.  The  wavelet  theory  proves  that  this  evo- 
lution across  scales  characterizes  the  local  Lipschitz  regularity  of  the  signal.  The  local  maxima  ot 
the  wavelet  transform  thus  not  only  detect  the  signal  sharper  variations  but  also  characterize  their 
"local  shape".  We  do  believe  that  this  complement  of  information  is  very  important  for  paitcm 
recognition.  An  "edge  detection"  should  not  only  find  the  signal  sharper  variation  points  but  also 
characterize  these  variations  in  order  to  distinguish  the  different  type  of  edges  and  discriminate 
textures.  We  describe  an  algorithm  to  compute  descriptors  of  the  edges  types  and  further  discuss 
their  applications. 

The  first  coherent  multiscale  edge  representation  was  developed  by  Marr  [24]  who  called  it 
"primal  sketch".  The  edges  of  the  primal  sketch  are  detected  from  the  zero-crossing  of  the  image 
convolved  with  the  Laplacian  of  a  scaled  Gaussian.  He  conjectured  that  such  a  representation  can 
be  complete  but  did  not  give  any  further  evidence.  The  mathematical  characterization  of  signals 
from  multiscale  zero-crossings  or  maxima  is  a  difficult  non-linear  mathematical  problem.  Wc 
review  some  important  results  and  explain  why  they  do  not  answer  Marr's  conjecture.  We  refor- 
malize  the  completeness  problem  within  the  wavelet  framework  and  show  that  it  leads  naturally 
to  an  algorithm  for  reconstructing  the  signal  from  its  multiscale  maxima.  This  algorithm  is 
described  and  implemented  both  for  one-dimensional  and  two-dimensional  signals.  It  provides 
for  the  first  time  an  exact  reconstruction  of  images  from  their  multiscale  edges  and  thus  verifies 
Marr's  conjecture. 

The  reconstruction  algorithm  shows  numerically  tiiat  the  whole  signal  information  can  be 
decomposed  in  multiscale  edges.  This  representation  has  potential  applications  in  signal  process- 
ing in  general  but  we  here  concentrate  on  image  processing  and  computer  vision.  Wc  discuss 
more  particularly  the  compact  coding  of  images.  The  wavelet  transform  maxima  not  only  detects 
edge  but  also  characterizes  their  local  shape.  This  opens  a  new  approach  to  the  regularization  of 
ill-defined  computer  vision  problems  without  making  any  smoothness  assumption  on  the  image. 
The  application  to  edge  detection  and  textures  discrimination  are  further  studied.  For  pattern 
recognition  purposes,  we  also  explain  how  to  integrate  coarse  to  fine  search  strategies  to  the 
wavelet  maxima  representation. 

We  introduce  the  most  important  mathematical  results  in  the  text  but  leave  the  proofs  in 
appendices.  The  transition  between  models  based  on  functions  of  continuous  variables  to  discrete 
computations  is  often  delicate  and  can  induce  many  errors.  We  tried  to  propcriy  build  this  transi- 
tion while  leaving  in  appendices  all  the  discrete  algorithms. 


Page  5 

Notation 

L^(R)  denotes  the  Hilbert  space  of  measurable,  square-integrable  one-dimensional  funciions 
fix).  For/(x)G  L  (R)  and  ^(x)  e  L  (R),  the  inner  product  of/ (x)  with  gU)  is  written: 

<g(.x),fix)>  =   J  six) fix)  dx. 

The  norm  (energy)  of/ (x)  e  L  (R)  is  given  by 

ll/lP  =  2l/(;c)|2  dx. 

We  denote  the  convolution  of  two  functions  fix)e  L  (R)  and  gix)  e  L'(R)  by 

+00 
/  *gix)  =    [  fiu)gix-u)  du. 

We  denote  by  H  (R)  the  Hilbert  space  of  differentiable  functions  in  the  sense  of  Sobolev.  The 
norm  of/  (;c )  g  H  (R)  is  given  by 

11/11^.   =  j^\fix)\^dx    +7|^i^|2^    . 

The  Fourier  transform  of  fix)  e  L^(R)  is  written /(co)  and  is  defined  by 

+00 
/(CO)  =  j  fix)e-''^  dx. 

L^(R^)  is  the  Hilbert  space  of  measurable,  square-integrable  two  dimensional  functions  fix,y). 
The  norm  of/ ix,y)e  L  (R^)  is  given  by: 

The  Fourier  transform  of/(x,>')  e  L^(R^)  is  wrinen  / (cOi  ,cOy )  and  is  defined  by 

/(CO,.a>,)  =  ^Jfix,y)e-'^'^-'^'^dx  dy  . 
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2.  Wavelet  Transform  Model  in  One  Dimension 

2.1.  Continuous  Wavelet  Transform 

We  briefly  introduce  the  wavelet  transform  and  describe  its  main  properties.  For  a  thorough 
presentation  the  reader  is  referred  to  a  general  review  [23]  and  an  advanced  functional  analysis 
book  of  Meyer  [27].  The  wavelet  model  has  first  been  fonmalized  by  Grossmann  and  Morlct 
[11].  Let  us  denote  by  \]isix)  the  dilation  of  \\iix)  by  a  factor  s: 

V.(j:)=y  V(j)  •  (1) 

The  wavelet  transform  of  a  function  fix)  at  the  scale  s  and  position  x  is  given  by  the  convolu- 
tion product: 

Wsf  (X)  =  f  *  xvsix)  .  (2) 

Let  y(a))  be  the  Fourier  transform  of  \\i{x).  Morlet  and  Grossmann  [11]  showed  that  if  ihc 
wavelet  \\i(x)  has  a  Fourier  transform  equal  to  zero  at  (o  =  0,  then  the  wavelet  transform  satisfies 
an  energy  conservation  equation  and  f{x)  can  be  reconstructed  from  its  wavelet  transform.  The 
wavelet  \^(x)  can  be  interpreted  as  the  impulse  response  of  a  band-pass  filter  and  the  wavciei 
transform  as  a  convolution  with  a  band-pass  filter  which  is  dilated.  If  the  scale  s  is  large, 
^sf(.x)  detects  the  lower  frequency  components  of  the  signal  f(x).  When  the  scale  s  decreases, 
the  support  of  \j/i(j:)  decreases  so  the  wavelet  transform  W^fix)  is  sensitive  to  firmer  details. 
The  scale  s  characterizes  the  size  and  the  regularity  of  the  signal  features  extracted  by  the 
wavelet  transform. 

22.  Dyadic  Wavelet  Transform 

The  wavelet  transform  depends  on  two  parameters  s  and  x  which  vary  continuously  over 
the  set  of  real  numbers.  For  practical  applications  these  parameters  must  be  discretized.  For  a 
particular  class  of  wavelets,  the  scale  parameter  can  be  sampled  along  the  dyadic  sequence 

2J  2.'  without  modifying  the  overall  properties  of  the  transform.  The  principles  of  such  a 
dyadic  scale  decomposition  was  studied  in  mathematics  by  Littlewood  &  Paley  in  the  1930's. 
The  wavelet  transform  at  the  scale  2-'   is  given  by: 

W^^fix)  =  f  *\V2,ix)  .  (3) 

At  each  scale  2-'.  the  function  W2,f(x)  is  continuous  since  it  is  equal  to  the  convolution  of  two 
functions  in  L  (R).  The  Fourier  transform  of  Wj,/^)  is 

Wf2,((i>)  =/((!))  vi;r(2V(o)  .  (4) 
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By  imposing  that 


S    |\)/(2^co)|2  =  1  . 


(5) 


we  insure  that  the  whole  frequency  axis  is  covered  by  a  dilation  of  \j/(a))  by  the  scales  factors 


21 


^g2-  ■^y  wavelet  satisfying  equation  (5)  is  called  a  dyadic  wavelet.  We  also  call  dyadic 


wavelet  transform  the  sequence  of  functions 


;6Z 


(6) 


We  denote  by  W  the  dyadic  wavelet  operator  defined  by  W/  = 


W^^fix) 


j^Z 


From  equations  (4)  and  (5)  and  by  applying  the  Parseval  theorem,  we  obtain  an  energy  con- 
servations equation: 

ll/l|2  =    I^I|W2,/(x)||2  .  (7) 

Let  \i/2;U)  =  '^2i^-x).  The  function  f  {x)  is  reconstructed  from  its  dyadic  wavelet  transform, 

f{x)  =     I^W2;/*V2;(^)   •  (8) 

This  equation  is  proved  by  computing  its  Fourier  transform  and  inserting  equations  (4)  and  (5). 


Let  V  be  the  space  of  the  dyadic  wavelet  transforms    VVjz/U) 


J7€Z' 


for  all  the  functions 


f(,x)e  L  (R).    Let  us  denote  by    I  (L  )    the  Hilbert  space  of  all  sequences  of  functions 


gjix) 


;eZ 


,  such  that 


gj{x)BL\R)   and       g    IU,U)l|2<+o= 


Equation  (7)  proves  that  V  is  a  sub-space  of  l^(L^).  We  denote  by  W  '   the  operator  from 
l^(L^)  to  L^(R)  defined  by: 


w- 


(9) 


The  reconstruction  formula  (8)  shows  that  the  restriction  of  W"^  to  the  wavelet  space  V  is  the 
inverse  of  the  dyadic  wavelet  transform  operator  W. 

Any  sequence  of  functions      gjix)    ^^e  F(L^)    is  not  a  priori  the  dyadic  wavelet 

transform  of  some  function /U)  g  L^(R).  Indeed,  if  there  exists  a  function  /U)  €  L"(R)  such 


that 


gjix) 


J€Z 


=  W/ ,  then  cleariy  we  should  have: 


W 


w 


-1 


gjix) 


jeZJ 


gjix) 


jeZ   ■ 


(10) 
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If  we  replace  the  operators  W  and  W~'  by  their  expression  given  in  equations  (3)  and  (9),  wc 
obtain; 

Vy  €  Z      ^f^g,*  Kijix)  =  gjix)  ,  with  (11) 

The  sequence     gjix)  is  a  dyadic  wavelet  transform  if  and  only  if  equations  (11)  hold. 

These  equations  are  called  reproducing  kernel  equations.  They  express  the  correlation  between 
the  functions  W2,f(x)  of  a  dyadic  wavelet  transform.  The  correlation  between  the  functions 
^2//^^^  ^""^  ^Tffi^)  at  two  different  scales  2-'  and  2'  can  intuitively  be  understood  by  look- 
ing at  their  Fourier  transform; 

W2,fi(0)  = /(co)v(2^co)    and    W.Jidi)  = /(co)  \j/(2'(o)  . 

The  redundancy  of  VVj/ZCco)  and  VVji/Cco)  depend  upon  the  overlap  of  the  functions  V|/(2^(o) 
and  \j/(2'a)).  The  energy  of  this  overlap  is  equal  to  the  energy  of  the  kernel  Kij(x).  It  is  max- 
imum for  I  =j  -  I  and  I  =j  +  I.  The  operator 

Pv  =  W  0  W-'  (12) 

is  a  projector  from  1  (L  )  on  the  V  space.  Indeed,  one  can  easily  prove  that  any  sequence  of 


functions 


gj{x)     ^2  €  1  (L  )  satisfies  Py  gj{x)     ^^  £  V  ,  and  we  saw  that  any  element  of 

V  is  invariant  under  the  action  of  this  operator.  The  projector  Py  is  important  for  the  purpose  of 
this  paper. 

2-3.  Wavelet  Transform  Amplitude  Across  Scales 

An  important  property  of  the  wavelet  transform  is  its  ability  to  characterize  the  local  regu- 
larity of  signals.  In  this  section,  we  explain  how  the  amplitude  of  the  wavelet  transform  relates  to 
the  singularities  appearing  in  the  signal.  The  local  regularity  of  a  function  at  a  point  xq,  can  be 
measured  with  the  Lipschitz  criteria.  Let  us  recall  the  definition  of  the  Lipschitz  regularity.  Let 
a  >  0  and  «  e  Z  be  such  that  n  <  a<«-t-I.  A  function  /(a:)  is  Lipschiu  a  regular  in  jtq  if 
and  only  if  there  exists  a  polynomial  Pni.x)  of  order  n  such  that  for  all  x  in  a  neighborhood  of 
xq,  we  have 

\f{x)-Pn{x)\    =    Oi\X-Xon    .  (13) 

Clearly, /(x)  is  then  n  times  continuously  differentiable  in  zq  and  the  polynomial  P„{x)  is 
the  first  n  terms  of  the  Taylor  series  of /(j:).  The  Lipschitz  regularity  a  thus  gives  an  informa- 
tion on  the  differentiability  of  the  signal  but  it  is  more  precise.  If /(a;)  is  n    times  continuously 
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differentiable  but  not  n+\  times,  the  coefficient  a  characterizes  the  singularity  of  the  n'^ 
derivative.  We  call  a  singularity  of  /(x),  any  point  xq  where  the  function  /U)  is  not  continu- 
ously differentiable.  In  a  singular  point,  the  Lipschitz  regularity  is  smaller  than  1  so  equation  (13) 
can  be  rewritten: 

\f{x)-f{xo)\  =  0(U-;cor).  with   a<l    .  (14) 

The  exponent  a  measures  the  "regularity"  of  the  singularity  in  ;co  and  characterizes  locally  the 
shape  of  the  function  around  the  singiilarity.  Let  Oq  <  0,  and  let  us  suppose  that  fix)  is 
Lipschitz  a  in  xq  if  and  only  if  a  <  Oo-  If  0  <  cxo  <  1,/(j:)  is  continuous  inxo  but  its  deriva- 
tive is  not  defined.  If  ao  =  0, /(a:)  is  bounded  but  discontinuous  in  j:o-  If  Oo  <0  ,/(x)  is  not 
bounded  in  xq.  In  order  to  understand  the  signification  of  the  Lipschiu  regularity  for  a  <  0  ,  one 
can  observe  that  if  a  function  is  Lipschitz  a  in  xq,  then  its  primitive  is  Lipschitz  a-i-1  at  the  same 
point.  The  primitive  of  a  Dirac  is  a  bounded  discontinuous  function  which  is  locally  Lipschitz  0. 
A  Dirac  in  xq  is  therefore  Lipschitz  -1.  The  following  theorem  relates  the  Lipschitz  regularity  of 
a  functions  to  its  wavelet  transform  amplitude.  We  suppose  that  the  wavelet  m^{x)  is  continu- 
ously differentiable  and  that   |  y(;c )\=0 (       ;) . 

Theorem  1 

Let  /(x)e  L^(R)./(x)  is  Lipschitz  a  in  all  points  of  an  open  interval  if  and  only  if  for  all  x  in 
this  interval 

\W^f{x)\  =  0(2«0  .  (1-^) 


The  proof  of  this  theorem  can  be  found  in  several  articles  [12, 15].  This  theorem  relates  the 
Lipschitz  regularity  of  a  function  to  the  asymptotic  decay  of  the  wavelet  transform  modulus  when 
the  scale  2'  decreases  to  zero.  This  characterization  is  valid  for  open  intervals  but  not  isolated 
points.  Jaffard  [15]  showed  that  the  regularity  of  a  function  in  a  particular  point  can  also  be 
derived  asymptotic  decay  of  the  wavelet  transform  in  the  neighborhood  of  this  point,  but  this 
results  is  slightly  more  technical.  Theorem  1  is  however  enough  for  estimating  the  Lipschitz 
regularity  of  isolated  singularities.  Let  us  suppose  that  /(a:)  has  a  Lipschitz  a  singularity  in  a 
point  xq  and  that  there  exists  an  open  interval  where  /(x)  is  continuously  differentiable  in  all 
points  different  from  xq.  In  this  interval  (including  xq),  f  (x)  is  Lipschitz  a  so  the  wavelet 
transform  satisfies  equation  (15).  The  coefficient  a  can  therefore  be  obtained  from  the  decay  of 
the  wavelet  transform  amplitude.  The  ability  of  the  wavelet  transform  to  characterize  the  type  of 
local  singularities  is  a  major  motivation  for  its  application  to  detect  the  signal  sharper  variations. 
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For  digital  processing  applications,  the  spatial  parameter  x  of  the  functions  V^'2^/  (x )  must 
also  be  discretized.  In  the  next  paragraph,  we  briefly  review  the  classical  linear  discretization  of 
a  wavelet  transform.  We  also  explain  the  translation  problems  which  appear  with  linear  discreti- 
zations. 

2.4.  Wavelet  Orthonormal  Bases  and  Translation  Invariance 

The  wavelet  transform  of  a  function  f(x)  at  the  scale  2-'  is  defined  in  equation  (3)  by  a 
convolution  product  However,  Wj^/Cx)  can  also  be  written  as  an  inner  product  in  L"(R): 

^/ 2/(^0)  =  /  *  V2/(^o)  =    j  f(x}\\i2,(x-xo)dx   .  thus 

W/2;Uo)   =    <f{x),\il2,(x-Xo)>    ■ 

The  wavelet  transform  of  fix)  at  the  scale  2^  at  ;io  is  therefore  equal  to  the  inner  product  of 
fix)  with  the  function  \\i2,ix-xo).  Meyer  [28]  and  Stromberg  [34]  have  shown  independently 
that  for  some  particular  choices  of  the  wavelet  \\iix),  the  family  of  functions 
\\f2jix-2Jn)  is    an   orthonormal    basis    of    L  (R)    .    Consequently,    the    samples 


,     ,  _,   provide  a  discrete,  uncorrelated  and  complete 


W:^fi2Jn)  =  <fix).yif^ix-2Jn)> 

representation  of  the  function  /  ix ).  We  showed  how  these  orthononmal  bases  model  the  concept 
of  multiresolution  and  we  described  some  applications  to  image  processing  [21].  At  first  glance 
one  could  think  that  such  a  representation  is  well  suited  for  pattern  analysis.  This  is  not  the  case 
because  of  its  non-invariance  with  translations.  Let  /(x)eL*'(R)  and  f-:ix)=fix-x)  be  a 
translation  of  fix).  Since  a  wavelet  transform  at  a  scale  2-'  is  given  by  a  convolution  product 
(equation  (3)),  it  is  clear  that  V^2'/t(x)  =  W2,fix-z)  .  However,  if  one  samples  each  of  theses 
functions  at  the  rale  2~J ,  the  values  obtained  may  be  totally  different  for  W^jfix)  and 
W2jf-:ix),  unless  'i  =  k2J  with  /:  e  Z  (see  fig.  1).  When  a  pattern  translates,  the  onhogonal 
wavelet  coefficients  are  therefore  modified.  It  is  thus  difficult  to  perform  any  pattern  recognition 
using  these  coefficients. 

A  possible  solution  to  the  translation  problem  is  to  oversample  each  function  W2,fix)  in 
order  to  obtain  a  quasi-translation  of  the  samples.  This  .solution  increases  considerably  the 
amount  of  data  needed  for  representing  the  signal  and  we  obtain  a  representation  that  is  very 
redundant.  The  representation  is  then  heavy  to  process  and  it  is  difficult  to  separate  die  properties 
of  the  signal  from  the  artifacts  due  to  the  inner  correlation  of  the  decomposition. 

In  order  to  obtain  a  translation  invariant  discretization  of  the  wavelet  transform,  one  can 
also  make  an  adaptive  sampling  of  each  function  Wjifix).  Qearly,  the  zero-crossings  and  the 
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"2./(x) 


W,JA^) 


Fig.  1:  This  drawing  shows  that  the  sampling  of  a  wavelet  transform  (given  by  the  arrows)  can  be 
very  different  after  translating  the  signal.  If  fxix)=f(x-x)  then  Wj^/tU)  =  Wj/U-t)  but 
sampling  does  not  translate  if  x  is  not  proportional  to  2> . 

local  extrema  of  the  functions  Wj^/Cj:)  do  translate  when/U)  translates.  The  zero-crossings 
and  local  extrema  of  H'2j/(;c)  are  well  defined  since  ^2// (^ )  is  continuous.  For  some  particu- 
lar wavelets  that  we  describe  in  the  next  section,  these  particular  points  locate  the  sharper  varia- 
tion points  of  the  signal.  The  two-dimensional  extension  of  this  procedure  is  equivalent  to  a  mul- 
tiscale  edge  detection  as  defined  in  computer  vision.  The  next  section  describes  the  propenies  of 
this  edge  detection  procedure. 


2J.  Wavelet  Extrema  and  Zero-Crossings 

The  local  extrema  and  zero-crossings  of  the  wavelet  transform  have  similar  properties. 
Their  interpretation  depends  upon  the  choice  of  the  wavelet  Mf{x ).  Let  us  call  a  smoothing  func- 
tion, the  impulse  response  of  a  low-pass  filter.  It  is  a  function  whose  Fourier  transform  has  an 
energy  concentrated  in  the  low-frequencies.  Let  e(z)  be  such  a  smoothing  function.  A  classical 
example  often  used  in  computer  vision  is  the  Gaussian.  Let  v/'C^)  and  vj/^U)  be  respectively  the 
first  and  second  order  derivative  of  Q{x): 


y(;c)  =  -^   and   y^Hx)  =  ^^ 


dx       —    -^  '"'         dx' 
The  wavelet  transforms  defined  with  respect  to  each  of  these  wavelets  are  given  by: 
Wjjfix)  =  f  *  \\ii(x)   and    Wifix)  =  f  *  vj/^U)  . 


(16) 


(17) 


Let  us  show  that  for  any  scale  2J  ,  the  zero  crossings  of  Wlf  (x )  correspond  to  the  local  extrema 
of  Wlf(x).  Let  Q2,ix)  =  jj-Q(,^)  be  the  dilation  of  Q{x)  by  a  factor  2J  .  As  a  consequence 
of  equations  (16),  one  can  rewrite  the  wavelet  transforms  W^fix)  and  Wlf{x): 

dQ. 


Wifix)  =f*{2>  -3^)ix)  =  ^'  -^(f*  92.)(^)    and 


(18) 
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Wlfix)  =  /  *  (22v  ^)(x)  =  22/  -^(f  *  Q^){x)  .  (19) 

The  wavelet  transforms  W^fix)  and  W^f{x)  are  proportional  respectively  the  first  and  second 
derivative  of /(a:)  smoothed  by  Bj^U)-  The  local  extrema  of  Wlfix)  thus  correspond  to  the 
zero-crossings  of  W^fix)  and  to  the  inflection  points  of  /  *  02;^)  •  The  extrema  detection  of 
the  wavelet  transform  built  from  \\i\x)  is  essentially  equivalent  to  a  Canny  edge  detection  [3] 
whereas  the  zero-crossing  detection  in  the  case  of  \\i\x)  can  he  viewed  as  a  Marr-Hildroih  [25  j 
edge  detector.  Kronland-Martinet  et.al.  [17]  have  developed  another  approach  for  delecting 
abrupt  changes  with  complex  wavelets.  This  approach  has  some  similarities  to  the  zero-crossing 
procedure  and  is  particularly  well  suited  for  complex  signals.  We  do  not  discuss  this  method  any 
further  but  rather  concentrate  on  the  comparison  between  zero-crossing  and  local  extrema. 

Detecting  zero-crossings  or  local  extrema  are  similar  procedures  but  the  local  extrema 
approach  has  some  important  advantages.  An  inflection  point  of  /  *  Qijix)  can  either  be  a  max- 
imum or  a  minimum  of  the  absolute  value  of  its  first  derivative.  The  maxima  of  the  absolute 
value  of  the  first  derivative  are  sharp  variation  points  of  /  *  Qi'^^^  whereas  the  minima 
correspond  to  slow  variations  (see  fig.  2).  These  two  types  of  inflection  points  can  be  dis- 
tinguished by  looking  whether  an  extremum  of  |W^/(x)|  is  a  maximum  or  a  minimum  but 
they  cannot  be  differentiated  from  the  zero-crossings  of  W^f  (x).  Since  we  want  to  detect  sharp 
variation  points,  we  only  record  the  maxima  of  |W^/(;c)|.  One  can  also  observe  that  zero- 
crossings  by  themselves  do  not  provide  a  stable  signal  representation.  Indeed,  a  slight  perturba- 
tion of  fix)  in  a  region  where  |  W^f  (x )  |  has  a  small  amplitude,  may  considerably  modi  fy  the 
zero-crossing  positions.  It  is  thus  necessary  to  stabilize  such  a  representation  with  a  complemen- 
tary information.  This  can  be  achieved  for  example  by  computing  the  wavelet  transform  energy 
between  two  consecutive  zero-crossings.  This  approach  has  already  been  studied  [20],  but  it  com- 
plicates slightly  the  representation.  When  detecting  a  local  maximum  of  \Wlf(x)\,  we  can 
record  both  its  position  xq  and  the  amplitude  Wlf(xo).  This  amplitude  stabilizes  the  represen- 
tation and  is  particularly  important  to  describe  the  signal  singularities.  Indeed,  we  saw  in  section 
2.3  that  the  evolution  of  the  wavelet  transform  amplitude  across  scales  characterizes  the  Lipschitz 
regularity  of  the  signal.  To  simplify  the  text,  in  the  rest  of  the  paper  we  say  that  we  detect  the 
maxima  of  the  wavelet  transform  although  we  really  mean  the  points  where  the  absolute  value  of 
the  wavelet  transform  is  maximum. 
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/(-x) 


/*6>2.(x) 


W-2^/(x) 


T^l./(x) 


Fig.  2:  r/ie  extrema  of  VV^/(a:)  a«d  f/ze  zero-crossings  of  W^f{x)  appear  at  the  inflection 
points  of  f  *  62;U).  We  only  record  the  points  of  abscissa  xq  and  X2  where  I  W;|,/(a:)|  is 
maximum  because  they  locate  the  sharper  variation  points  of  fix)  smoothed  at  the  scale  2' , 
The  local  minima  of  |W^/(x)|  in  x\  also  corresponds  to  a  zero-crossing  to  Wl,f(x)  but  it 
locates  a  slow  variation  point. 


3.  Completeness  and  Stability 

A  fundamental  issue  is  to  understand  whether  the  local  maxima  or  the  zero -crossings  derinc 
a  complete  and  stable  representation  of  the  original  signal.  In  section  3.1,  we  review  some  results 
obtained  on  similar  problems.  Most  previous  works  have  been  done  with  zero-crossings.  We  saw 
in  section  2.5  that  zero-crossings  and  local  extrema  provide  the  same  type  of  information  depend- 
ing upon  the  choice  of  the  wavelet  These  studies  can  thus  be  extended  from  zero-crossings  lo 
extrema. 


3.1.  Previous  Mathematical  Results 

The  most  classical  result  concerning  the  characterization  of  a  signal  from  its  zero-crossings 
is  due  to  Logan  [18].  We  describe  in  some  detail  Logan's  theorem  because  it  provides  a  good 
understanding  of  the  mathematical  issues.  Let  f{x)e  L^(R)  and  let  us  suppose  thai  its  Fourier 
transform  has  a  support  included  in  one  octave  intervals.  Logan  theorem  [18]  proves  that  if 
fix)  does  not  share  any  zero-crossings  with  its  Hilbert  transform,  then  it  is  uniquely  character- 
ized by  its  zero-crossings.  Let  us  give  an  intuitive  justification  of  this  result.  Wc  know  that  there 
exists     0)0     such    that    the    Fourier    of    fix)     has    a    support    included    in    the    intervals 
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[-2coo  .  -coq]  \j  [ooo  .  2(Oo]-  The  Nyquist  theorem  proves  that  such  a  signal  is  characterized  by  a 

uniform  sampling  at  the  rate   — ^  .  One  can  also  prove  that  this  signal  changes  sign  approxima- 

tively  as  frequently  as  the  function  sin(coox:)  .  The  number  of  zero-crossings  is  therefore  of  the 
same  order  than  the  number  of  values  needed  to  characterize  the  signal  with  a  uniform  sampling. 
Of  course,  the  zero-crossing  problem  is  different  since  zero-crossings  are  not  uniformly  distri- 
buted but  one  can  see  that  qualitatively  the  same  amount  of  information  is  available.  To  prove 
this  theorem,  Logan  makes  an  analytical  extension  of  the  signal  and  uses  standard  properties  on 
zeros  of  analytical  functions.  The  zero-crossing  characterization  as  explained  by  Logan  is  not 
stable:  "the  problem  of  actually  recovering  (the  signal)  from  its  sign  changes  appears  lo  be  very 
difficult  and  impractical". 

Let  us  now  explain  how  Logan's  theorem  can  be  integrated  in  the  wavelet  model.  Let  \\i(x ) 
be  the  function  equal  to  the  impulse  response  of  a  perfect  band  pass  filter  of  one  octave  ifs 
Fourier  transform  is  given  by: 

-     ^      J  1  if  7t<  |co|  <27: 
V(w)  =  I  0  otherwise  (20) 

The  function  \)/(x)  clearly  satisfies  equation  (5)  and  is  therefore  a  dyadic  wavelet.  Let 
fix)e  L^(R);  the  Fourier  transform  of  W2,fix)  is  given  by  W 2,/ (oi)  =  f  ((o)  \\i(2J (a)  .  The 
support  of  W2if(.(i))  is  thus  included  in  the  one  octave  intervals 
[-2~J'*'^K  ,-2~Jk]\^[2~Jk  ,2~J'^^k]  .  From  Logan's  theorem  we  derive  that  each  function 
^2if(x)     is    characterized   by    its    zero-crossings.    Since    we    can    reconstruct    f(x)     from 

W^fix) 
functions 


,  the  original  function  f{x)  is  also  characterized  by  the  zero-crossings  of  all  the 

je  Z 


W^^fix) 


.  This  characterization  is  however  not  stable  as  previously  explained. 


Although  Logan's  theorem  is  an  important  result,  we  want  now  to  emphasize  the  reason 
why  it  cannot  be  used  for  the  type  of  wavelets  we  are  interested  in.  We  need  a  wavelet  equal  to 
the  second  derivative  of  some  smoothing  function  so  that  zero-crossings  indicate  the  position  ol 

the  signal  shaper  variation  points.  If  \\i(,x)  =         y^  ,  then  its  Fourier  transform  can  be  wnttcn 

Y((o)  = -<o^  6(co) .  Since  6(z)  is  a  smoothing  function,  it  satisfies  0(0)  ?tO  so  \)/((o)  has  a  zero 
of  order  two  at  co  =  0.  Similariy,  one  can  show  that  a  wavelet  is  the  n'''  order  derivative  of  some 
smoothing  function  only  if  its  Fourier  transform  \|/(g))  has  a  zero  of  order  n  at  co  =  0.  The 
Logan  wavelet  \\i(x)  given  in  (20)  cannot  be  written  as  a  finite  order  derivative  of  some  smooth- 
ing function  since  its  Fourier  transform  has  an  infinite  order  zero  in  (0  =  0.  Hence,  the  zero- 
crossings  of  the  wavelet  transform  W2jf(x)  can  not  be  interpreted  as  any  particulariy  interesting 
features  of  fix).   In  fact,  there  are  too  much  zero-crossings  since   W2,f(x)   changes  sign  in 
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almost  all  intervals  of  length  V ,  independently  of  f  {x).  Logan  as  well  as  other  researchers  who 
extended  this  result,  use  the  band-limited  properties  of  the  signal  for  computing  its  analytical 
extension.  Ail  these  proofs  do  not  provide  any  stability  result  since  they  are  based  on  non-stable 
characterization  of  analytical  functions  [5.33,38].  The  reader  is  referred  to  a  review  by  Hummel 
and  Moniot  for  more  details  [14]. 

Many  suidies  have  also  described  the  properties  of  zero-crossings  of  functions  convolved 
with  the  Laplacian  of  a  Gaussian.  This  convolution  is  equivalent  to  a  wavelet  transform  buili 
from  a  wavelet  v|/(j:)  equal  to  the  Laplacian  of  a  Gaussian.  Such  a  wavelet  transform  can  be  inter- 
preted as  the  result  of  a  heat  diffusion  process  [16].  Indeed,  the  Gaussian  is  the  Green  function  of 
the  heat  diffusion  equation.  Let  t  =s^  be  the  diffusion  time,  one  can  show  that  the  wavelet 
transform  W^f  {x)  built  with  the  Laplacian  of  a  Gaussian  satisfies  the  heat  differential  equation 

aw./U)    ^   d'WJ{x)  ^21) 

The  wavelet  transform  V/^fix)  is  therefore  equal  to  a  heat  distribution  after  a  diffusion  time 
r  =5^  with  an  initial  heat  distribution  at  r  =0  equal  to  A/U)  (the  Laplacian  is  taken  in  the 
sense  of  distributions).  Using  the  properties  of  the  heat  diffusion  differential  equations  several 
authors  have  proved  interesting  properties  of  the  propagation  of  zero-crossings  across  scales 
[37, 13, 16].  Hummel  &  Moniot  as  well  as  Yuille  &  Poggio  have  also  proved  that  the  position  of 
the  zero-crossings  of  W^fix)  give  a  complete  characterization  of  the  function  f{x)[\'i].  These 
proofs  are  based  on  deconvolution  arguments  which  are  highly  instable  and  they  do  not  take  a  full 
advantage  of  the  information  given  by  the  zero-crossings  at  aU  scales.  They  can  therefore  not  be 
used  for  the  reconstruction  of /(x)  from  the  zero-crossings  of  its  wavelet  transform.  The  dif- 
ferential equation  (21)  gives  the  evolutionary  properties  of  WJ{x)  when  the  .scale  s  and  ihc 
abscissa  x  vary.  It  expresses  the  redundancy  of  the  functions  WJ{x)  at  different  scales.  This 
differential  equation  fonnalism  has  several  inconveniences.  It  is  only  valid  for  a  wavelet  equal  lo 
the  Laplacian  of  a  Gaussian.  This  is  quite  restrictive  since  the  wavelet  transform  built  with  the 
Laplacian  of  a  Gaussian  cannot  be  implemented  with  a  fast  pyramidal  algorithm.  The  Laplacian 
pyramid  of  Burt  [2]  approximates  the  Gaussian  kernel  with  a  function  of  compact  support  no 
more  than  twice  continuously  differentiable.  It  is  not  clear  what  properties  derived  from  the  heal 
diffusion  equation  are  valid  for  such  a  kernel.  Also,  the  heat  fonnalism  cannot  be  extended  to  a 
discrete  sequence  of  scales.   In  all  applications,  the  scale  parameter  varies  on  sparse  discrete 


sequences  such  as  the  dyadic  sequence 


2' 


y.   One  cannot  use  the  heat  partial  differential 


equation  to  express  the  redundancy  specifically  at  these  scales. 

The  wavelet  transform  reproducing  kernel  is  an  alternative  approach  to  fonnalize  the  redun- 
dancy of  the  functions   WJ{x).  Grossmann  [11]  showed  that  the  wavelet  transfomi  satisfies  a 
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reproducing  kernel  equation 


w, 


■fix)  =    I  j^Kis.s'^,y)W,'f(y)dyds' 


(22) 


where  the  kernel  K(s,s'  jcjc')  depends  upon  the  wavelet  \\i(x).  When  the  wavelet  is  equal  to 
the  Laplacian  of  a  Gaussian,  this  equation  is  an  integral  equation  associated  with  the  heat  paniai 
differential  equation  (21).  However,  the  reproducing  kernel  equation  is  valid  for  any  wavelet 
\|/(x).  The  energy  of  the  kernel  K{sy  jc,y)  measures  the  redundancy  between  the  function 
Wsfix)  and  Ws-fix)  at  two  different  scale  s  and  s' .  The  reproducing  kernel  formalism  can 
be  extended  to  a  discrete  sequence  of  scales.  We  saw  in  section  2.2  that  for  dyadic  scales,  the 
functions  W2,f{x)  also  satisfy  a  reproducing  kernel  equation  given  by  equations  (11).  In  the 
next  section,  we  reformalize  the  completeness  problem  by  using  the  wavelet  transform  reproduc- 
ing kernel.  As  explained  in  section  2.5,  we  detect  the  local  maxima  of  the  wavelet  transform 
rather  than  the  zero-crossings  and  therefore  use  a  wavelet  equal  to  a  first  order  derivative  of  a 
smoothing  function. 


3.2.  Signal  Reconstruction  from  the  Wavelet  Maxima 

In  section  2.5,  we  differendated  the  maxima  of  \W2jf(x)\  from  its  minima  since  wc  arc 
interested  in  the  sharper  variation  points  of/  (x ).  We  first  study  the  reconstruction  of  fix)  from 
all  extrema  (maxima  and  minima)  of  \W2,fix)\  and  then  consider  the  case  where  only  the 
maxima  are  kept.  We  reformalize  the  completeness  problem  within  the  wavelet  framework  and 
then  derive  an  algorithm  to  perform  the  reconstruction. 


Let  fix)e  L^(R)  and 


Wz^fix) 


jeZ 


be  its  dyadic  wavelet  transform.  Since  f(x)  can 


W.Jix) 


^,^  given 


be  recovered  from  its  dyadic  wavelet  transform,  we  first  try  to  reconstruct 

the  local  extrema  of  each  funcnon  Wj^/U),  j  e  Z.   Qeariy,  for  any  scale   2^ ,  there  exists  an 

infinite  number  of  functions  gj(x)  which  have  the  same  local  extrema  as  Wt^/Cj:).  However, 


any  such  sequence  of  funcdons 


gjix) 


j;ez 


is  not  necessarily  the  dyadic  wavelet  iransfonn  ol 


some  function  in  L  (R).  Indeed,  we  saw  in  section  2.2  Uiat  for  being  a  dyadic  wavelet  transfonm, 
it  must  satisfy  the  reproducing  kernel  conditions  (11).  We  thus  have  two  types  of  information  for 


reconstructing  the  functions 


W2,fix) 


;eZ' 


We  know  the  positions  and  amplitudes  of  the  local 


extrema  of  each  function  VV2,/(x)  and  we  want  to  reconstruct  a  sequence  of  functions  which 
satisfies  the  inner  redundancy  given  by  the  reproducing  kernel  equations  (11).  Let  us  recall  from 


secdon  2.2  that    l^(L^)    is  the  space  of  all  sequence  of  functions      gjix) 


"^  \\gj(.x)\\'^  <-h=°  .  The  space  of  all  dyadic  wavelet  transforms 


W2,fix) 


je/. 


„    such  thai 

yeZ 

is  denoted   V 
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and  is  a  sub-space  of  1^(L  ).  We  suppose  that  the  wavelet  \\i(x)  is  differentiable  in  the  sense  of 
Sobolev:  \\i(x)  e  H'(R)-  At  each  scale  2J ,  we  then  have  W2,f(x)  €  H  (R).  In  order  to  express 
the  conditions  given  by  the  extrema  of  the  wavelet  transform  of  /(.x),  we  dcHnc  the  set  F  of  all 


sequences 


gjix) 


jeZ 


in  P(L')  such  that  for  all  scales  2-',  gj{x)e  H'(R)  and  gjix)  has 

the  same  extrema  than  W2jf{x).  We  prove  in  appendix  5  that  the  set   V  is  convex.  The  local 
extrema  representation  is  complete  if  and  only  there  exists  no  dyadic  wavelet  transform  different 


from 


Wj^fix) 


V  is  reduced  to  one  element: 


that  has  the  same  local  extrema.  In  other  words,  the  intersection  of  V  with 


n 


V  = 


W:^fix) 


jeZ 


(23) 


In  order  to  verify  numerically  this  assertion,  we  describe  an  algorithm  that  reconstructs  the  inter- 
sections of  r  with  V. 

A  classical  technic  for  recovering  the  intersection  of  a  convex  set  with  a  linear  space  is  to 
iterate  on  alternative  projections  on  the  convex  and  the  linear  space.  Youla  and  Webb  [36j  wrote 
a  review  of  the  mathematical  properties  of  these  algorithms.  Let  us  suppose  that  the  convex  F 


is  closed  within  an  appropriate  Hilbert  space.  For  any 


gjix) 


jeZ 


in  this  Hilbert  space,  we  can 


define  a  projection    Pp   on    F   that  transforms      gj{x) 


jsZ 


into  the  sequence  of  functions 


hjix) 


jeZ 


6  F  that  is  the  closest  to 


gj(x) 


jjeZ- 


Since  F  is  convex,  the  projection   Pp  is 


nonexpansive.  Let  Py  be  the  orthogonal  projection  on  the  space  V  and  P  =  Pp  o  Py  be  the 
composition  of  Pp  and  Py.  Qeariy  any  element  at  the  intersection  of  F  and  V  is  a  fixed 
point  of  P.  To  compute  such  a  fixed  point,  we  iterate  on  the  operator  P.  Let  P*"'  be  the  compo- 
sition n  times  of  the  operator  P.  Since  Pp  is  a  nonexpansive  projection  on  a  closed  convex  and 
Py    is  an  orthogonal  projection,  one  can  prove    [36]  that  for  any  sequence  of  functions 


8jM 


;ez 


,  when    n    tends  to    -h»  P^") 


gjix) 


jeZ 


converges  weaidy  to  an  element  in 


r^  V.  This  iterative  algorithm  reconstructs  the  dyadic  wavelet  transform  of  f(x),  whatever 


the  initial  sequence 


gjix) 


jeZ 


,  if  and  only  if  the  intersection  of  F  and  V  is  reduced  to  one 


element.  The  reconstruction  algorithm  is  illustrated  in  fig.  3. 
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Fig.  3:  The  reconstruction  of  the  wavelet  transform  of  f  {x)  from  its  extrema  (or  maxima)  con- 
sists in  iterating  on  the  operator  P.  This  operator  is  a  composition  of  the  projection  Pp  on  the 
set  T  which  expresses  the  constraints  on  the  extrema  {or  maxima],  and  the  projection  Pv  on 
V  which  is  the  space  of  all  the  dyadic  wavelet  transforms.  The  goal  is  to  converge  to  the  inter- 
section of  T  and  V. 


In  order  to  define  the  projector  Pf,  we  must  insure  that  F  is  closed.  The  convex  F  is  not 
closed  in  the  Hilbert  space  l'(L')  but  it  can  be  closed,  as  suggested  by  Daubcchics  [fi],  within  a 
Hilbert  space  derived  from  the  Sobolev  space  H'(R).  Appendix  5  characterizes  more  precisely 
this  Hilbert  space  and  the  corresponding  operator  Pp.  In  section  2.2.  we  saw  that  the  reproducing 
kernel  of  a  dyadic  wavelet  transform  defines  a  projection  operator  Pv  on  V.  Since  the  projec- 
tion operator  Py  is  given  by  a  reproducing  kernel,  one  can  prove  that  it  defines  an  onhogonal 
projection  on  V.  An  iteration  on  these  two  operators  Pp  and  Pv  thus  converges  to  an  element 
inFnV. 

If  instead  of  recording  the  extrema  of  Wj/U),  we  only  keep  the  maxima  of  \W2,f{x)\, 
we  can  build  a  similar  reconstruction  algorithm.  The  only  difference  is  the  definition  of  the  sci 
F.  In  this  case,  F  is  the  set  of  all  sequences  gj(x)\  ^^  in  1  (L*)  such  that  for  all  scales  V, 
gj(.x)€  H'(R)  and  gj{x)  has  the  same  maxzma  than  W2,f{x).  Since  we  loose  the  constraint 
on  the  minima  of  |  W-^f  {x)\,  we  prove  in  appendix  5  that  F  is  not  convex  any  more.  Wc  can 
still  define  an  operator   Pp   that  transforms  any  sequence  of  functions    Igjix)     ^.^    into  the 

However,  since  F  is 


sequence  of  functions    hj{x)\        e  F  that  is  the  closest  to  \gj{x)     ^.^. 

not  convex.  Pp  is  not  a  nonexpansive  operator.   Hence,  the  iteration  on  the  operator   P  is  not 

guaranteed  to  converge. 
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In  the  following  sections,  we  develop  the  tools  necessary  for  the  numerical  implementation 
of  this  algorithm.  The  algorithm  is  implemented  in  the  case  where  we  only  keep  the  maxima  of 
the  wavelet  transfoirn.  Although  we  have  no  proof  of  convergence,  the  numerical  experiments 
described  in  section  4.3  show  that  the  algorithm  converges  quickly  and  that  it  does  reconstruct 
the  original  wavelet  transform. 

4.  Discrete  One-Dimensional  Model 

A  proper  implementation  of  the  maxima  representation  and  the  reconstmction  algorithm 
raises  several  important  questions.  The  input  signal  is  generally  measured  with  a  finite  resolution 
which  imposes  a  firmer  scale  when  computing  the  wavelet  transform.  In  practice,  the  scale 
parameter  must  also  vary  on  a  finite  range.  Section  4.1  explains  how  to  interpret  mathematically  a 
dyadic  wavelet  transform  on  a  finite  range  of  scales.  In  all  previous  sections,  our  model  was 
based  on  functions  of  a  continuous  parameter  x.  We  discretize  the  abscissa  x  and  describe  effi- 
cient algorithms  for  computing  a  discrete  wavelet  transform  and  its  inverse.  Section  4.2  analyzes 
the  potential  errors  introduced  by  the  discretization  when  detecting  Uie  position  and  amplitude  of 
the  wavelet  transform  maxima.  The  numerical  implementation  of  the  algorithm  which  recon- 
structs the  signal  from  its  wavelet  transfomn  maxima  is  described  in  section  4.3  and  we  study  the 
numerical  results.  Section  4.4,  gives  a  procedure  to  compute  descriptors  of  the  signal  shaper  vari- 
ations from  the  amplitude  of  the  wavelet  maxima  at  different  scales. 


4.1.  Discrete  Dyadic  Wavelet  Transform 

In  practice,  we  cannot  compute  the  wavelet  transform  at  all  scales  2-'  for  j  varying  from 
-oo  to  -K».  We  are  limited  by  a  finite  larger  scale  and  a  non-zero  finner  scale.  Let  us  suppose  for 
normalization  purposes  that  Uie  finner  scale  is  equal  to  1  and  that  2-^  is  the  largest  scale.  Let 
/(x)g  L^(R).  We  first  show  that  between  the  scales  1  and  2\  tiie  wavelet  transform 
^2;/(^)  </  can  be  interpreted  as  tiie  details  available  when  smoothing  f{x)  at  the  scale  1 
but  which  have  disappeared  when  smootiiing  f{x)  at  tiie  larger  scale  2-'.  Let  us  introduce  a 
function  (t)(x)  whose  Fourier  transform  is  given  by: 


l<^((0)|2  =  ^  IV(2^a))|2   .  (25) 

Since  the  wavelet  v(x)  satisfies    X    I  vi/(2-' (o)  |  ^  =  1 ,  one  can  derive  tiiat  lim|$(CL))|  =  1  .  The 

energy  of  the  Fourier  transform   $(co)   is  concentrated  in  tiie  low  frequencies  so   <S)(x)   is  a 
smoothing  function.  Let  us  define  tiie  smootiiing  operator  ^-^  by: 

52./U)  =  /  *  <{>2/(^)  .  witii    '1)2.  =  ^  <f(^)  ■  (^6) 
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The  larger  the  scale  2-',  the  more  details  of  f(x)  are  removed  by  the  smoothing  operator  5 2 


Let  us  prove  that  the  dyadic  wavelet  transform 


W:^f(x) 


\ij^ 


between  the  scales  1  and  2^  pro- 


vide the  details  available  in  S\f(,x)  but  not  in  S2jf(x).  The  Fourier  transform  of  S]f{x)  , 
Sjjfix)  and  W2jf(x)  are  respectively  given  by 

S'lfm  =  $(co)/((o)    .   Syfim)  =  (^(2^(0) /(CO)  and                           (27) 

H/2,/((o)  =  v(2^co)/((o)  .  (28) 

Equation  (25)  yields 

l$(co)|2  =   |j|vi}(2^(o)|2  +   |$(2^(0)|2  .  (29) 

Using  Parseval's  theorem,  we  derive  from  equations  (27),  (28)  and  (29)  the  following  energy  con- 
servation equation: 


l|5i/(x)||2=|^  \\W-^f{x)\\^+\\Syf{x)\\- 


(30) 


This  equation  proves  that  the  higher  frequencies  of  Sif{x)  which  have  disappeared  in  Syf  {x ) 


can  be  recovered  from  the  dyadic  wavelet  transform 


Wyfix) 


\<j<j 


between  the  scales  1   and 


i<j<j 


are  called  \he  finite-scale  wavelet  transform  of 


2^  .  The  functions  ISyfix),    W^f  {x ) 

Sif(x).  In  practice,  the  signal  we  process  is  given  by  a  discrete  sequence  of  values.  The  follow- 
ing lemma  proves  that  any  discrete  signal  of  fmite  energy  can  be  interpreted  as  the  uniforan  sam- 
pling of  some  function  smoothed  at  the  scale  1. 


Lemma  1 


Let  D  = 


d„        „  be  a  discrete  signal  of  finite  energy:     Y    |c/„r  <  -h=o.  Let  us  suppose  thai 
the  Fourier  transform  <f(G))  satisfies: 

3Ci>03C2>0,  with   V(0€  R     ,     Ci  <    J    | $(a)+2n tt) | ^  <  d  ■ 

There  exists  a  function  f(x)e  L  (R)  (not  unique)  such  that 

V/iG  Z  .  SJ{n)  =  d,    .  (31) 


pies 


The  proof  of  this  lemma  is  in  appendix  1.   The  discrete  signal    D  can  thus  be  rcwriticn 
.  For  a  particular  class  of  wavelets  \\i{x)  described  in  appendix  2,  ihc  sam- 

enable  us  to  compute  a  uniform  sampling  of  the  finite  scale  wavelet 


D=    Sifin) 


Sxfin) 


transform  of  Sif(x) 


Page  21 


Syfin) 


neZ  ' 


W^f(n+w) 


neZ 


\<.j<J 


(32) 


The  sampling  shift  w  depends  upon  the  wavelet  V|/(a:).  Let  us  denote 


^if  =     W^f(n+w) 


neZ 


and   Sif  =    S^f(^n) 


neZ 


(33) 


The  sequence  of  discrete  signals 


SU  , 


Wif 


\<.j^ 


is  called  a  discrete  dyadic  wavelet 


transform  of  the  signal  D  = 


Sxfin) 


„.  We  denote  by  W*^   the  discrete  wavelet  transform 

operator  which  associates  to  a  signal  D  the  discrete  wavelet  transform  previously  defined. 
Appendix  3  describes  a  fast  algorithm  for  implementing  this  operator.  If  the  signal  has  A'  non 
zero-samples,  the  complexity  of  this  algorithm  is  O  {N  log{N)).  It  is  based  on  a  cascade  of  con- 
volutions with  two  discrete  filters  H  and  G .  Appendix  3  also  describes  the  implementation  of 
the  discrete  inverse  wavelet  transform  W"'-^  which  reconstructs  the  signal  D  from  its  discrete 
dyadic  wavelet  transform.  The  reconstruction  algorithm  has  also  a  complexity  of  O  (N  log  (N)). 

Fig.  4(a)  is  the  graph  of  a  cubic  spline  wavelet  described  in  apjjendix  2.  The  signal  in  fig. 
5(a)  is  an  image  scan-line  of  256  samples  and  fig.  5(b)  is  its  discrete  wavelet  transform  computed 
with  the  cubic  spline  wavelet.  The  curves  in  both  figures  are  linear  interpolations  between  the 
samples  of  each  discrete  signals.  The  curve  at  the  bottom  of  fig.  5(b)  is  the  coarse  signal  S^J . 
Since  the  wavelet  used  is  the  first  derivative  of  a  smoothing  function,  the  maxima  of  the  wavelet 
transform  indicate  the  points  of  sharper  variation  at  each  scale. 


Fig.  4:  (a):  Graph  of  a  cubic  spline  wavelet  ^(x).  Its  computation  is  described  in  appendix  2.  It 
decreases  exponentially  at  infinity. 

(b):  Graph  of  the  function  i^{x)  corresponding  to  the  wavelet  \\i{x)  plotted  in  fig.  4(a).  This  func- 
tion is  also  a  cubic  spline  but  has  a  compact  support  of  size  4. 
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(c) 
Fig.  5:  Each  of  these  graphs  are  linear  interpolation  Oetween  the  samples  of  discrete  signals. 
(a):  Image  scan-line  of  256  samples. 

(b):  Dyadic  wavelet  transform  of  signal  5(a)  decomposed  on  5  scales. 

(c):  Maxima  representation  of  signal  5(a).  At  each  scale,  the  Diracs  indicate  the  position  and 
amplitude  of  a  local  maximum  of  the  wavelet  transform  given  in  5(b). 
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A2.  Maxima  Detection  of  a  Discrete  Wavelet  Transform 

Lemma  1  proves  that  any  discrete  signal  is  equal  to  the  unifomi  sampling  of  some  function 


Sifix)=f  *  ^(x) .  We  can  thus  define  the  discrete  wavelet  transform  <  S^,f  , 


wif 


\<i<j 


which  is  a  uniform  sampling  of  the  finite  scale  wavelet  transform  of  f  {x).  At  each  scale  2i , 

1  <  y  <y ,  we  need  to  compute  the  maxima  of  V/2,f  {x)  from  its  uniform  sampling  W^/ .  Let  us 

suppose  for  simplification  that  the  sampling  shift  w   is  zero.   We  call  discrete  maxima  of  the 


discrete       signal 


wi  = 


W^fin) 


such        that 


,       any       sample        VVi,/(n) 

Wvf{n)\  >  \W2,f{n-\)\  and  \W2,f{n)\  >  \W2,f{n+\)\  .  We  make  no  errors  by  estimat- 
ing the  maxima  of  W2,f(x)  from  the  discrete  maxima  of  W^f  if  and  only  if  each  maxima  of 
^2//^^^  h^  ^  integer  abscissa.  The  following  lemma  gives  a  sufficient  condition  to  satisfy 
this  maxima  property. 


Lemma  2 

Let   us   suppose   that   the   function    ^{x)    is   such   that   there   exists    p{x)e  L'(R)    such 

S ip(x)  =  p  *  ^(x)  satisfies 


5ip(x)  =  Sip(-x)      Vj:  G  R 
Sipix)  =  0    if    |jc|  >1 

^-^^iP^<0     if   0<;c<l 
dx 

S  ipix)  +  S  ipi\-x)  =  I     ifO<x<l 


For  any  discrete  signal  D,  there  exists  f(x)e  L  (R)  such  that  D  - 
scale  2-'  >  1  the  maxima  of  Wj^/Cx)  have  integer  abscissa. 


(34) 


SJin) 


leZ 


and  for  all 


The  proof  of  lemma  2  is  written  in  appendix  4.  The  function  5ip(x)  is  a  non-oscillatory 
interpolation  of  a  discrete  Dirac: 

^     ^  ^      J  1  if  «=0 
Sip(n)  =  |o  if  /I  ^0  ■ 

We  can  express  the  function  Sifix)  of  lemma  2  from  the  sample  values  of  D  and  S  ,p(A  ).  in 
order  to  find  such  a  function  p(x)  the  support  of  <\)(x)  must  at  most  have  a  size  2.  Indeed,  the 
support  of  Sip(x)  is  always  larger  than  the  one  of  <t)(x)  and  condition  (34)  implies  that  it  has  a 
size  2.  A  particular  example  of  function  ^{x)  where  there  exists  p{x)  which  satisfies  the  condi- 
tions (34)  of  lemma  2  is 


_  I   1  if  0  <  X  <  1 
*'*(-^)  ~  I  0  otherwise 


(35) 
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The  corresponding  dyadic  wavelet  is  the  Haar  wavelet: 


\if{x)  = 


1  if  0<x  <  ^ 

-1  if  ^  <  Jt  <  1  .  (36) 

0  otherwise 


We  can  verify  that  if  p(x)  =  ^(x+\)  then  Sip{x)  satisfies  (34)  and 

c  «.  ^      J  1-UI   if  Ul  <1  ,„, 

S:Pix)  =|o  if  |z|>l  (37) 

The  wavelet  maxima  representation  of  a  signal  D  is  defined  by  the  discrete  maxima  of  its 
discrete  wavelet  transform  W^f  for  each  scale  \<2'<2-',  plus  the  coarse  signal  Syf  ■  Fig. 
5(c)  gives  the  maxima  representation  of  the  discrete  wavelet  transform  shown  in  fig.  5(b).  Each 
Dirac  indicates  the  position  and  amplitude  of  a  maxima.  As  it  can  be  observed,  these  ma.xima 
indicate  at  different  scales  the  position  of  the  signal  sharper  variations.  We  show  in  appendix  4 
that  the  functions  ^(x)  whose  support  is  larger  than  2  do  not  satisfy  lemma  2  and  in  this  case  the 
discrete  maxima  detection  introduces  errors.  We  explain  that  the  errors  can  be  interpreted  as  the 
addition  of  spurious  high  frequencies  in  the  signal  that  we  decompose.  For  example,  the  cubic 
spline  function  ^{x )  shown  in  fig.  4(b)  has  a  support  larger  than  2. 

43.  Numerical  Reconstruction  from  the  Wavelet  Transform  Maxima 

The  algorithm  that  reconstructs  the  original  signal  from  its  local  maxima  represemaiion,  is 
based  on  two  projection  operators.  The  first  one  is  the  projection  Py  on  the  space  of  all  dyadic 
wavelet  transforms.  To  any  sequence  of  functions,  it  associates  the  dyadic  wavelet  transform  of 
some  function  f  (x)  e  L^(R).  We  saw  in  equation  (12)  that  this  operator  can  be  decomposed  into 

Pv  =  W  0  W'    . 

where  W  and  W~'  are  respectively  the  wavelet  and  inverse  wavelet  transform.  Within  a 
discrete  framework,  this  operator  is  redefined  by 

P^  =  W  0  W-i-^  . 

where  W  and  W"'-^  are  respectively  the  discrete  wavelet  transform  and  the  inverse  discrete 
wavelet  transform.  Since  both  W'^  and  W"'-"^  are  implemented  with  fast  algorithms  of  com- 
plexity 0{N  log{N)),  the  numerical  complexity  for  implementing  P^  is  also  0{N  log{N)). 
The  other  projection  operator  involved  in  the  reconstruction  is  the  non-linear  projection  on  the  set 
r.  Appendix  5  describes  the  discrete  implementation  of  this  operator  that  we  denote  Pf^ .  The 
implementation  of  Pf^  has  a  complexity  0(N  log{N}).  Let  P^  =Pf^o  P^  .  The  reconstruction 
algorithm  iterates  on  the  operator  P'  to  reconstruct  the  intersection  of  V  and  F.  We  begin  the 
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iteration  from  an  initial  sequence  of  discrete  signals  built  with  a  spline  interpolation  between  the 
maxima  of  the  wavelet  transform.  We  then  iterate  on  the  operator  V^  until  it  converges  to  a 
fixed  point.  We  apply  the  inveree  wavelet  transform  operator  W"'-^  in  order  to  compute  the  sig- 
nal corresponding  to  the  reconstructed  dyadic  wavelet  transform. 


s^^J 


so 


lOO  150 

(a) 


200 


25  O 


Fig.  6:  (a):  Reconstruction  of  the  dyadic  wavelet  transform  from  the  maxima  representation  given 
in  fig.  5(c).  This  reconstruction  was  obtained  after  15  iterations  on  the  operator  P^ . 
(b):  Reconstruction  of  the  signal  by  applying  the  inverse  wavelet  operator  W''-^   on  the  recon- 
structed wavelet  transform  of  fig.  6(a).  The  quality  of  the  reconstruction  can  be  appreciated  by 
comparing  this  graph  with  fig.  5(a).  The  noise  to  signal  ratio  of  the  reconstruction  is  2  10"^. 
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Fig.  6(a)  shows  the  reconstruction  of  the  dyadic  wavelet  transform  from  the  local  maxima 
representation  given  in  Fig.  5(c),  with  15  iterations  on  the  projector  P^.  Fig.  6(b)  is  the  recon- 
struction of  the  original  signal  by  applying  the  inverse  wavelet  transform  operator  on  fig.  6(a). 
The  error  signal  is  computed  by  subtracting  the  reconstructed  signal  from  the  original  one.  Wc 
compute  the  ratio  of  the  signal  error  energy  to  the  energy  of  the  signal  minus  its  mean  and  the 
square  root  is  called  noise  to  signal  ratio.  We  suppress  the  mean  when  computing  the  noise  to 
signal  ratio  because  it  is  not  reconstructed  from  the  maxima  but  it  is  kept  in  the  coarse  signal 
Syf .  The  noise  to  signal  ratio  of  the  reconstruction  given  in  fig.  6(b)  is  2  10"^.  Fig.  7(a,b)  give 
the  noise  to  signal  ratio  for  the  reconstruction  of  the  discrete  wavelet  transform  signals 
W^f  ,  as  a  function  of  the  number  of  iterations  on  the  operator  P^.  For  each  scale  2-* , 

the  noise  to  signal  rauo  does  not  converge  to  zero  but  to  a  minimum  asymptotic  error.  The  noise 
to  signal  ratio  of  the  remaining  error  decreases  approximatively  by  a  constant  factor  when  the 
scale  increases  by  2.  It  proves  that  the  error  is  concentrated  in  the  high  frequencies.  This  error 
can  be  explained  by  the  spurious  high  frequencies  introduced  by  the  maxima  detection  procedure. 
Indeed,  these  numerical  experiments  are  performed  with  the  cubic  spline  wavelet  of  fig.  4(a)  and 
we  saw  in  the  previous  section  that  in  this  case,  the  maxima  obtained  from  a  discrete  wavelet 
transform  are  not  exactly  equal  to  the  maxima  of  a  finite  scale  wavelet  transform.  Appendix  4 
shows  that  the  corresponding  error  can  be  viewed  as  spurious  high  frequencies  that  are  added  to 
the  signal.  Fig.  7(c,d)  plot  the  noise  to  signal  ratio  of  the  signal  obtained  by  applying  the  inverse 
wavelet  transform  operator  on  the  reconstructed  discrete  wavelet  transforms.  After  10  iterations, 
the  noise  to  signal  ratio  already  reaches  2.4  10"^.  However,  because  of  the  remaining  high  fre- 
quency errors,  the  noise  to  signal  ratio  does  not  converge  to  zero  but  to  a  minimum  value  approx- 
imatively equal  to  8  \0~^. 

The  reconstruction  algorithm  has  been  tested  on  a  large  collection  of  functions:  Diracs, 
sinusoidal  waves,  edges,  singularities,  fractal  curves,  image  scan  lines...  For  aO  theses  tests  the 
ratio  of  the  error  to  signal  energy  has  approximatively  the  same  behavior  than  in  fig.  7.  In  particu- 
lar, the  error  is  monotonously  decreasing  as  if  the  operator  P^  was  nonexpansive.  This  result  was 
not  predicted  by  the  analysis  in  section  3.2  since  the  set  V  is  not  convex  when  we  keep  only  the 
maxima  of  the  wavelet  transform  instead  of  all  extrema.  The  high  frequency  reconstruction  error 
introduced  with  the  cubic  spline  wavelet  is  small  enough  for  many  applications.  However,  wc 
can  suppress  the  remaining  high  frequency  errors  by  choosing  another  wavelet  such  as  the  Haar 
wavelet  so  that  the  corresponding  function  ^(x)  satisfies  lemma  2.  In  this  case,  we  prove  in 
appendix  4  that  the  maxima  detection  is  exact.  Fig.  8(a)  gives  the  noise  to  signal  ratio  when 
reconstructing  the  discrete  wavelet  transform  of  signal  5(a)  from  the  wavelet  maxima  obtained 
with  a  Haar  wavelet.  With  a  Haar  wavelet,  the  reconstruction  error  decrease  up  to  the  limit  of  the 
computation  precision.  After  4000  iterations,  the  signal  to  noise  ratio  is  smaller  than  10"^  at  all 
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Fig.  7:  (a),(b):  The  ordinate  gives  the  noise  to  signal  ratio  when  reconstructing  the  discrete 
wavelet  transform  signals  from  the  wavelet  maxima  of  fig.  5(c),  on  a  logarithmic  scale.  Each 
curve  labeled  by  a  number  (k)  gives  the  noise  to  signal  ratio  for  the  reconstruction  of  Wjj .  The 
abscissa  is  the  number  of  iterations  on  the  operator  P^ .  These  experiments  are  done  with  the  cu- 
bic spline  wavelet  of  fig.  4(a).  The  noise  to  signal  ratio  does  not  converge  to  zero  because  of  the 
high  frequency  errors  introduced  by  the  maxima  detection. 

(c),(d):  The  ordinate  gives  the  noise  to  signal  ratio  of  the  signal  obtained  from  the  reconstructed 
wavelet  transform,  on  a  logarithmic  scale.  The  abscissa  is  the  number  of  iterations  on  the  opera- 
tor P^. 


scales.  After  100  iterations,  the  decay  rate  of  the  noise  to  signal  ratio  is  slow  but  approximaiively 
constant.  For  a  given  number  of  iterations,  the  signal  to  noise  ratio  decreases  approximatively  by 
a  constant  factor  when  the  scale  increases  by  2.  This  indicates  that  the  remaining  error  arc  con- 
centrated in  the  high  frequencies.  The  reconstruction  of  the  higher  frequencies  is  slower  because 
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of  the  structure  of  the  reproducing  kernel  for  discrete  signals.  We  saw  in  section  2.2  that  the 
reproducing  kernel  expresses  the  redundancy  of  the  wavelet  transform  at  different  scales  and  that 
the  redundancy  between  Wj^/Cx)  and  WjJix)  is  maximum  if  /  =)  -  1  or  /  =7  +  1.  For  a 
discrete  wavelet  transform,  at  the  scale  2',  the  reproducing  kernel  expresses  the  redundancy 
between  W^^f  and  Wj,f  but  there  is  no  constraint  with  the  finner  scale  2^  since  we  can  not 
compute  iL  The  reproducing  kernel  thus  imposes  less  constraint  on  the  high  frequencies  which 
explains  why  their  reconstruction  is  slower.  Fig.  8(b)  gives  the  noise  to  signal  ratio  for  the  recon- 
struction of  the  signal  itself.  After  4000  iterations,  the  noise  to  signal  ratio  is  smaller  than 
2  10^.  The  noise  to  signal  ratio  had  the  same  behavior  for  all  other  experiments  perfoirned  with 
the  Haar  wavelet  This  seems  to  indicate  that  for  this  wavelet,  the  maxima  representation  is  com- 
plete and  stable.  It  is  impxjrtant  to  observe  again  that  the  noise  to  signal  ratio  is  monotonically 
decreasing  as  if  the  operator  P^  was  nonexpansive.  In  practice,  the  number  of  iterations  thai  is 
need  depends  upon  the  reconstruction  precision  that  is  needed.  The  value  2  10"^  is  well  below 
the  measurement  precision  of  most  signals. 
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Fig.  8:  (a):  The  ordinate  gives  the  noise  to  signal  ratio  when  reconstructing  the  discrete  wavelet 
transform  signals  from  the  wavelet  transform  maxima  of  signal  5(a),  computed  with  a  Haar 
wavelet.  Each  curve  labeled  by  a  number  (k)  gives  the  noise  to  signal  ratio  for  the  reconstruction 
of  W^,f .  The  abscissa  is  the  number  of  iterations  on  the  operator  P'^ .  The  noise  to  signal  ratio 
decreases  up  to  the  limit  of  the  computation  precisions. 

(b):  The  ordinate  gives  the  noise  to  signal  ratio  of  the  signal  obtained  by  applying  the  inverse 
wavelet  operator  on  the  reconstructed  wavelet  transform,  computed  with  a  Haar  wavelet.  The 
abscissa  is  the  ruunber  of  iterations  on  the  operator  P^ . 


All  the  nimierical  results  obtained  so  far.  seem  to  indicate  that  the  maxima  of  a  dyadic 
wavelet  transform   do   provide   a  complete   and  stable   representation.   The   high-frequency 
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reconstruction  errors,  when  using  a  cubic  spline  wavelet,  seem  to  be  due  to  the  maxima  detection 
errors  since  it  disappears  when  using  a  Haar  wavelet  Such  an  error  clearly  does  not  exist  when 
stating  the  problem  in  continuous  variables.  We  therefore  conjecture  that  for  a  large  class  of 


dyadic  wavelets,  the  maxima  of 


W:^nx) 


provide  a  complete  and  stable  represcniaiion  of 


/  (;c ).  The  class  of  wavelet  for  which  this  is  true  remains  to  be  defined.  We  want  to  stress  ihai  this 
is  only  a  conjecture  based  on  numerical  results  but  no  proof  is  given  in  the  paper.  The  recon- 
struction algorithm  previously  described  is  however  sufficient  for  the  applications. 

4.4.  Behavior  of  Maxima  Amplitude  Across  Scales 

The  previous  section  showed  experimentally  that  we  can  reconstruct  the  original  signal 
from  the  maxima  of  the  wavelet  transform  and  thus  that  this  representation  is  complete  and 
stable.  Let  us  now  explain  how  to  use  directly  the  maxima  representation  for  characterizing  the 
different  type  of  signal  variations. 

Theorem  1  proves  that  we  obtain  the  Lipschitz  regularity  of  a  signal  from  the  behavior  of  its 
wavelet  transform  amplitude  when  the  scale  decreases.    In  order  to  apply  this  theorem,  ihc 

wavelet  "^{x)  must  be  differentiable  and  satisfy    I v/(x)|  =  0{       j).  The  Haar  wavelet  cannoi 

be  used  since  it  is  discontinuous  but  we  can  choose  the  cubic  spline  wavelet  given  in  fig.  4(a). 
Let  us  first  analyze  the  behavior  of  the  wavelet  transform  maxima  when  the  signal  has  an  isolated 
singularity.  By  isolated  singularities,  we  mean  any  singular  point  xq  where  f  {x)  is  Lipschiu 
a  but  not  Lipschitz  oo  for  some  Oo  >  a,  whereas  f  {x)  is  Lipschitz  Oo  in  all  points  diffcrcni 
from  xq  in  a  neighborhood  of  j:o-  If  Oo  >  l./(^)  is  continuously  differentiable  in  the  neighbor- 
hood of  Xq  but  if  0(0  <  1 ,  f  (x)  might  be  singular  in  this  neighborhood  although  it  is  less  singular 
than  in  xq.  A  typical  example  is  a  step  edge  in  xq  (a  =  0)  with  a  texture  on  both  side  of  the  edge 
which  is  Lipschitz  Oo.  with  0(o  >  0.  Let  us  prove  that  the  singularity  of  f  {x)  in  xq  creates  local 
maxima  of  the  wavelet  transform  'W2,f{x)  for  scales  2'  small  enough,  in  the  neighborhood  of  xq. 
For  any  point  .xi  *a:o  in  a  neighborhood  of  xo,  we  can  find  an  open  interval  around  x  \  where  /  {x ) 
is  Lipschitz  ccq  in  all  points  of  this  interval.  Theorem  1  proves  that  in  this  interval  we  have 
I  ^2if  (^ )  I  =0  {2°^ ).  On  the  contrary,  in  any  open  interval  which  contains  xq,  since  /  (x )  is  not 
Lipschitz  Oo  in  xq,  theorem  1  proves  that  the  wavelet  u-ansform  does  not  satisfy  this  decay  pro- 
perty. One  can  thus  derive  that  for  some  scales  21  small  enough,  the  wavelet  transform  W2,f{x) 
has  at  least  one  local  maxima  in  the  neighborhood  of  xq.  When  the  scale  V  decreases  to  zero, 
the  positions  of  these  local  maxima  converge  to  xq.  Since  f  {x)  is  Lipschitz  a  in  a  neighbor- 
hood of  xo,  theorem  1  proves  that  the  amplitude  of  the  wavelet  transform  is  bounded  by  0(2"'' ) 
so  the  amplitude  of  its  maxima  is  also  bounded  by  O  (2"-J ). 
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Often  the  signal  singularities  are  blurred  due  to  some  diffusion  process.  We  thus  gci 
smoothed  "singularities"  and  it  is  important  to  estimate  this  smoothing  factor.  For  example,  the 
shadows  in  an  image  do  not  produce  sharp  discontinuities  of  the  image  intensity  but  relatively 
smooth  variations  because  of  the  diffraction  effect.  In  this  case,  the  image  intensity  is  locally 
equal  to  the  convolution  of  a  step  edge  with  a  smoothing  Icemel.  The  width  (scale)  of  the  smooth- 
ing kernel  depends  upon  the  amount  of  diffraction.  Let  us  prove  that  this  smoothing  component 
can  be  estimated  from  the  behavior  of  the  wavelet  transform  amplitude  across  scales.  Let  g  {x ) 
be  a  function  which  is  Lipschitz  a  in  a  given  neighboitiood.  We  suppose  gix)  is  locally 
smoothed  at  the  scale  ^o-  To  simplify  the  analysis,  we  suppose  that  the  smoothing  kernel  Q(x )  is 
the  primitive  of  the  wavelet  \\i(x).  Let  f(x)  =  g  *  QsJ.^)-  We  saw  in  equation  (18)  that  the 
wavelet  transform  of  f{x)  can  be  written: 

W:^f{x)  =  2J  ^(f  *  Q2,)ix)  =  2J  -^(g*  e,„  *  e.Jix) . 

Let  us  suppose  that  for  any  pair  of  scales  (so  ,2J ),  the  function  Q{x )  satisfies 

Qv  *64^)  =  ^s,ix)  .  (38) 

where  the  scale  ,$1  dep»ends  upon  {sq  ,  2'):  s\  =  (3(2-' Jo)-  Gaussians  for  example  satisfy  this 
equation  with  P(2-'  ^q)  =  \2'^J  +  io^  •  The  wavelet  transform  W2,f{x)  can  be  rewritten 

W^fix)  =  21  -^(^*  e,,)(;c)  =  H-  W,^g{x)  .  (39) 

In  theorem  1,  we  characterized  the  behavior  of  the  wavelet  transform  at  the  scale  1'  from  the 
Lipschitz  regularity  of  the  signal.  This  theorem  can  easily  be  extended  for  all  positive  scales. 
Since  g{x)is  locally  Lipschitz  a  in  an  open  interval,  we  can  prove  that 

\W,^g{x)\   =  0{sn  ■  (40) 

Equation  (40)  yields 

\W^f{x)\   =  0{2'  sr^)     with     sx  =  ^{2'^q)    .  (41) 

One  can  easily  show  that  for  any  smoothing  function  Q(x )  which  satisfies  (38),  we  have 


5,  =  (i(2^Jo)  =  ' 


2'   if  21  »  So 

So  if  so»2i    ■  ^"^^' 


We  can  thus  distinguish  two  extreme  domains.  If  2'  »  sq,  then  s\~2'  so 

\W^f{x)\  =  0{2'^)  .  (43) 

In  this  range  of  scale,  the  wavelet  transform  is  not  sensitive  to  the  smoothing  at  the  scale  so.  It 
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behaves  as  if  the  function  is  Lipschitz  a  regular.  If  2'  «  so,  then  s\~  sq  so 

\W^f{x)\  =  0(50"-' 2^)  =  OiU)  .  (44) 

In  this  range  of  scale,  the  wavelet  transform  W2,f(.x)  decays  like  2-'.  This  indicates  that  the 
wavelet  transform  is  not  sensitive  any  more  to  the  singularities  that  might  exist  in  g  {x ),  because 
of  the  smoothing  effect.  When  the  scale  2-'  is  of  the  order  of  the  smoothing  scale  sq,  then  both 
the  singularities  of  g(x)  and  the  smoothing  influences  the  wavelet  transform  amplitude. 

We  explained  in  section  4.1  that  in  practice  we  have  an  approximation  at  finite  resolution  of 
f(x)  and  can  only  compute  the  wavelet  transform  for  scales  larger  than  1  (after  renormaliza- 
tion).  Equation  (41)  gives  the  behavior  of  the  wavelet  amplitude  across  scales  when  a  singularity 
Lipschitz  a  is  smoothed  at  the  scale  sq-  We  can  estimate  the  parameters  a  and  ^o  by  assuming 
that  there  exists  a  constant  K  such  that  the  amplitude  aj  of  the  maxima  at  the  firmer  scales  2-' 
satisfy  exactly: 

aj=K2Jsr^    with   ii  =  (3(5 0.2^  .  (45) 

In  this  case,  by  analogy  with  equation  (41),  we  say  that  in  xq  the  signal  approximated  at  the  reso- 
lution 1  behaves  like  a  singularity  Lipschitz  a  smoothed  at  the  scale  ^o-  For  a  given  a  and  sq, 
the  constant  K  is  proportional  to  the  amplitude  of  the  signal  change.  For  a  step  edge  a  =  0  and 
K  is  proportional  to  the  height  of  the  edge.  In  order  to  compute  the  three  constant  K  ,  a  and  ^o. 
we  must  at  least  know  the  evolution  of  the  maxima  amplitude  across  three  scales  for  a  given 
singularity.  One  can  also  consider  more  scales  and  compute  the  three  constants  such  that 
^  2-'  sp-'  (s\  =  p(5o,2-') )  approximates  the  maxima  amplitude  at  each  scale,  with  a  minimum 
error.  Let  us  suppose  that  we  choose  compute  K  ,  a  and  sq  from  the  maxima  amplitude  at  /  dif- 
ferent scales  (/  >  3).  Such  a  computation  is  meaningful  only  if  the  the  maxima  detected  from  ilic 
scales  2'  to  2'  do  correspond  to  a  unique  singularity  that  we  want  to  characterize.  For  exam- 
ple, the  wavelet  transform  shown  in  fig  5(b)  has  a  unique  maxima  between  the  abscissa  60  and 
80  at  the  scale  2^.  However,  this  maxima  does  not  correspond  to  a  unique  isolated  singularity  but 
to  aU  the  irregularities  of  the  signal  5(a)  in  the  same  interval.  At  finner  scales,  we  can  indeed  ver- 
ify that  this  maxima  breaks  into  several  maxima  of  comparable  amplitude  which  converge 
towards  the  different  sub-variations  of  the  signal.  On  the  contrary,  in  the  abscissa  160  of  the 
same  figure,  the  signal  has  a  sharp  singularity  which  dominates  in  a  large  neighborhood.  The 
corresponding  maxima  between  the  scales  2'*  and  2'  do  not  break  into  several  maxima  of  com- 
parable amplitude.  Whether  a  singularity  is  isolated  or  not  at  a  given  scale  is  one  of  the  criteria  to 
distinguish  texmre  points  and  edge  points.  Typically,  the  signal  irregularities  of  fig.  5(a)  between 
the  abscissa  60  and  80  are  considered  as  a  "texture"  whereas  the  sharp  change  in  the  abscissa  160 
is  considered  as  an  "edge".  The  differentiation  between  textures  and  edges  and  the  discrimination 
of  textures  is  one  of  the  most  difficult  problem  of  low-level  image  processing.  We  further  discuss 
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this  issue  in  section  7.3. 


In  the  following,  we  concentrate  on  the  characterization  of  isolated  "edges".  We  study  the 
signal  singularities  which  create  extrema  on  more  that  3  scales  and  suppose  that  the  amplitude 
of  these  extrema  only  depends  upon  the  local  behavior  of  a  unique  singularity  and  not  of  several 
signal  variations.  In  order  to  compute  the  value  of  the  constants  K ,  a  and  so,  we  must  relate 
the  maxima  which  appear  at  different  scales  and  are  created  by  the  same  signal  change.  Lei 
(Xj,aj)  be  a  maximum  at  the  scale  2J  whose  position  is  Xj  and  amplitude  aj .  We  need  to  find 
the  maxima  at  the  scale  2-'"'  which  is  created  by  the  same  signal  variation.  Among  all  the  max- 
ima at  the  scale  2J~^  having  the  same  sign  than  Uj ,  we  select  the  ones  whose  abscissa  arc  closer 
to  Xj  than  any  other  maxima  at  the  scale  2-'.  Among  all  these,  we  select  the  maxima  whose 
amplitude  absolute  value  is  the  largest.  If  the  maxima  ixj,aj)  does  correspond  to  a  unique  signal 
change  and  not  to  an  agglomeration  of  several  ones,  the  algorithm  makes  the  right  correspon- 
dence. Otherwise  this  correspondence  is  not  meaningful  anyway  since  each  maxima  is  not  pro- 
duce by  one  but  several  signal  changes.  Let  us  suppose  that  we  related  three  maxima  appearing 
at  the  scales  2',  2^  and  2^  which  correspond  to  the  same  signal  change  and  let  us  compute  the 
constants  K ,  a  and  sq  of  this  signal  variation.  Let  ai,  aj  and  a^  be  respectively  the  amplitude  of 
these  maxima  at  the  scale  2',  2^  and  2^.  From  equation  (45),  we  derive  that  the  three  constants  K , 
a  and  ^o  of  this  signal  change  must  satisfy: 

ai=K  2^  P(2i,io)"-*  ,  a2  =  K2^  P(22,5o)"-'  and  03  =  K  2^  p(2l:fo)""'  ■  (46) 

These  three  equations  are  sufficient  to  compute  K ,  a  and  sq-  For  the  purpose  of  these  compula- 
tions, we  chose  a  wavelet  yU)  whose  primitive  Q(x)  is  similar  to  a  Gaussian  and  satisfies  with  a 
good  approximation: 


02.  *  6.0  =  0.,    w'tli    ^1  =  P(2^5o)  =  V22v  +S6    .  (47) 

This  wavelet  is  a  quadratic  spline  and  is  described  at  the  end  of  appendix  2.  Let  F{x)  be  the 
function  given  by 

Let  F~^(x)  be  the  inverse  of  F(x).  This  inverse  has  no  analytical  expression  but  is  easily  com- 
puted numerically.  With  a  few  algebraic  manipulation  on  equations  (46)  and  (47),  we  prove  that 

^o-'^     ^bgia2)-logia3)  +  logi2)^    '  ^^^^ 

a  =  2 -^^^llz^iip>±^  ^  1    and    K  =  ^  (2^ ^ sS)"^  .  (50) 

tog  (2^  +  50^)-  log  (2"  +  Sq^)  2 

We  applied  these  formula  to  estimate  the  Lipschitz  regularity  a  and  the  smoothing  scale  io  of 
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shaip  variation  points  from  the  behavior  of  their  wavelet  maxima  amplitude  across  scales.  Fig. 
9(a)  is  the  linear  interpolation  of  a  discrete  signal  Sff  which  is  buUt  from  a  function  /  (x ) 
which  has  four  singularities.  The  first  one  is  a  step  edge  (a  =  0)  whereas  the  second  one  is  a  step 
edge  smoothed  at  the  scale  .Jo  =  5.8.  The  third  singularity  is  Lipschitz  a  =  0.25  with  no  smooth- 
ing where  as  the  fourth  one  has  the  same  Lipschitz  regularity  but  is  smoothed  at  the  scale 
.so  =  5.8.  Fig.  9(b)  is  the  discrete  wavelet  transform  at  the  scales  2',  2^  and  2^.  Clearly  the 
behavior  of  the  maxima  at  different  scales  is  different  for  each  of  these  singularities.  The  table  1 
gives  the  results  of  numerical  experiments  for  the  estimation  of  the  Lipschitz  regularity  a  and 
the  smoothing  scale  sq  from  the  wavelet  maxima  at  the  scales  2',  2^  and  2^  We  take  four  dif- 
ferent singularities  whose  Lipschitz  regularity  a  is  respectively  equal  to  0  ,  0.25  ,  0.5  and  0.75. 
Each  of  these  singularities  are  smoothed  with  a  smoothing  scale  sq  respectively  equal  to 
0 , 2  , 3.5  and  5.8  .  For  each  resulting  signal  we  estimate  the  Lipschitz  regularity  a  and  the 
smoothing  scale  sq  from  the  maxima  of  the  wavelet  transform  by  applying  equations  (49)  and 
(50).  The  results  of  these  computations  are  reported  in  table  1.  The  value  of  the  Lipschitz  regu- 
larity a  is  recovered  with  a  good  precision  but  the  error  increases  when  the  smoothing  scale  sq 
increases.  The  error  on  the  estimation  of  sq  also  increases  with  sq.  These  errors  are  due  to  the 
fact  that  the  primitive  Qix)  of  the  wavelet  used  in  these  experiments,  satisfies  equation  (38)  only 
approximatively  and  the  approximation  error  increases  with  ^o-  One  can  take  into  account  this 
error  by  adding  other  terms  to  equation  (38)  and  get  more  precise  estimate  of  5o  and  a.  How- 
ever, these  results  already  show  that  the  estimation  of  the  parameters  a  and  so  is  precise  enough 
to  differentiate  many  types  of  singularities  from  the  behavior  of  the  wavelet  transform  maxima 
across  scales. 
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a. 

^0 

a 

^0 

0.00 

0.00 

0.00 

0.08 

0.00 

2.00 

0.02 

2.00 

0.00 

3.50 

0.07 

3.21 

0.00 

5.80 

0.14 

4.90 

0.25 

0.00 

0.25 

0.05 

0.25 

2.00 

0.26 

2.00 

0.25 

3.50 

0.30 

3.23 

0.25 

5.80 

0.35 

4.95 

0.50 

0.00 

0.50 

0.04 

0.50 

2.00 

0.51 

2.00 

0.50 

3.50 

0.53 

3.25 

0.50 

5.80 

0.56 

5.01 

0.75 

0.00 

0.75 

0.02 

0.75 

2.00 

0.75 

2.00 

0.75 

3.50 

0.76 

3.27 

0.75 

5.80 

0.78 

5.07 

Table  1:  The  first  two  columns  give  the  Lipschitz  regularity  a  and  the  smoothing  scale  sq  of 
several  singularities.  The  two  right  columns  are  the  estimate  of  the  Lipschitz  regularity  a  and 
the  smoothing  scale  sq'  obtained  by  applying  equations  (49)  and  (50)  on  the  amplitude  of  the 
wavelet  maxima,  created  by  each  singularity.  The  estimation  errors  increase  with  sq. 


stf 


(a) 


A 


-A. 


w'f 


(b) 


Fig.  9:  (a):  The  signal  Sif  is  built  from  a  function  f{x)  which  has  four  singularities  respec- 
tively characterized  by  {a  =  0  ,so  =  0)  ,  (a  =  O.Jo  =  5.8)  ,  (a  =  0.25 .5o  =  0 )  and 
(a  =  0.25,  Jo  =  5.8). 

(b):  The  behavior  of  the  maxima  of  W^,f  ,  V/iif  and  W^,f  depend  upon  the  Lipschitz  regularity 
a  and  the  smoothing  scale  sq. 
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5.  Wavelet  Maxima  Representation  of  Images 


5.1.  Dyadic  Wavelet  Transform  of  Two-Dimensional  Signals 

Let  us  now  extend  the  wavelet  local  maxima  model  to  two-dimensional  signals.  The  dyadic 
wavelet  transfomi  is  extended  in  L^(R2)  by  introducing  several  wavelets  ^''(x,y)  €  L'(R-)  , 
I  <k  <K .  Each  wavelet  ^'^{x,y)  can  have  a  specific  orientation  tunning.  Let  us  denote  by 

The  wavelet  transfomi  of  a  function  f(,x,y)e  L^CR^)  at  a  scale  2J ,  in  a  point  (x,y)  and  for  an 
orientation  k  is  defined  by 

Wif{x,y)=f  *^i{x,y)  .  (51) 

Let  ^P*  (©I  ,03, )  be  the  Fourier  transform  of  ^*  ix,y).  The  Fourier  transform  of  w!^f(x ,  v )  is 

Wif  ((o, .cOy )  =  f\(0, ,(0y )  ^* (U (0, ,2J(^)  .  (52) 

In  order  to  insure  that  we  can  cover  the  whole  frequency  plane  by  dilating  4'*(a)i,cOy)  by  the 


scales  factors 


2J 


g  2'  *£  impose  that 


f    g    |4>*(2^co,.2^o)y)|2  =  1 


(53) 


In  fig.  10,  each  wavelet  ^'^(x,y),  \  <  k  <K  ,  has  a  specific  orientation  tuning  and  the  suppon  of 
their  Fourier  transform  cover  together  approximatively  an  annulus.  In  two  dimensions,  wc  also 
denote  by  W  the  dyadic  wavelet  operator  which  associates  to  any  function  f{x,y)e  L  (R-^)  iis 
dyadic  wavelet  transform     Wl^f  (x  ,y ) 


jeZ.lSASAT' 


Fig.  10:  In  this  example,  the  frequency  supports  of  T*((0, ,C0y)  1  <^  <6  cover  an  annulus. 
Each  wavelet  has  a  specific  orientation  selectivity.  The  dilation  of  this  annulus  by  factors  2' 
should  cover  the  whole  frequency  plane  to  satisfy  the  condition  (53). 
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The  properties  of  a  two-dimensional  wavelet  transform  are  essentially  the  same  as  in  one 
dimension.  These  properties  are  proved  from  equations  (52)  and  (53)  and  Parseval's  theorem.  The 
energy  of  a  function  f  {x,y)  is  equal  to  the  sum  of  the  energies  of  the  functions  W^/  (x  ,>■ )  for 
all  the  orientations  k  and  scales  2' : 


(54) 


Let    'i'^{x,y)  =  ^l(-x,-y).   The  function /(x.y)    is  reconstmcted  from  its  dyadic  wavelet 
transform  by 

K 


f(x,y)  =  T  i:  wy  *^i{x,y)  . 


(55) 


We  denote  by  W~'  the  inverse  wavelet  transform  operator  which  reconstructs  f  {x,y)  from  its 
dyadic  wavelet  transform.  As  in  the  one-dimensional  case,  the  correlation  between  the  dyadic 
wavelet  transform  functions  can  be  expressed  by  reproducing  kernel  equations.  Any  sequence  of 


functions 


gjKx,y) 


jeZl^iK 


is  not  a  priori  the  two-dimensional  dyadic  wavelet  transform  of 


some  function  f(x,y)e  L^(R^).   There  exists  a  function  /(;c,>')e  L*(R2)   such  that  for  all 


jeZ\<k<J< 


IS  mvan- 


scales  2-*   and  orientations  k,  g/(a:,>')  =  W^/(j:,>'),  if  and  only  if     gjKx,y) 

ant  imder  the  action  successively  of  the  inverse  wavelet  transform  operator  followed  by  the 

wavelet  transform  operator: 


W 


w- 


g}ix,y) 


]€Z.\<Jc<Ji 


gJKx  ,y ) 


jeZ.  \<Jc<J< 

J    ■' 


Similarly  to  the  one-dimensional  case,  by  replacing  the  operator  W  and  W"'  by  their  analytical 
expression,  we  obtain  reproducing  kernel  equations.  Let  us  denote  by  V  the  vector  space  of  all 
two-dimensional  dyadic  wavelet  transform  of  functions  in  L  (R^).  The  operator 

Pv  =  W  0  W-'  (56) 

is  a  projector  on  V.  Since  this  projector  is  given  by  a  reproducing  kernel  equation,  it  is  onhogo- 
nal  with  respea  to  an  L^(R2)  norm  applied  on  the  sequence  of  functions     gf{x,y)       ■i\<k<K- 

In  two  dimensions,  the  regularity  of  a  function  fix.y)e  L  (R^)  can  also  be  derived  from 
the  evolution  of  the  wavelet  transform  amplitude  across  scales.  For  a<  1,  a  function  f(x.y)  is 
Lipschitz  a  in  (,xo,yo)  if  and  only  if 

\f(x,y)-f(xo,yo)\  =  0(U-xor+  \y-yon  ■  (57) 

If  a  function  is  not  continuously  differentiable  in  a  point  then  its  Lipschitz  regularity  in  this  poini 
must  be  smaller  or  equal  to  1.  A  simple  extension  of  theorem  1  proves  that  f{x,y)  is  Lipschii/ 
a  in  an  open  set  of  R^  if  and  only  if  for  all  points  {x  ,y )  in  this  set 
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\Wif(x,y)\   =  OH'")     for  \<k<K  .  (58) 

The  Lipschitz  regularity  can  therefore  be  obtained  from  the  asymptotic  behavior  of  the  wavelet 
transform  amplitude  when  the  scale  2J  tends  to  zero. 

For  digital  computations,  we  need  to  discretize  the  spatial  coordinates  (x,y).  A  possible 
approach  is  to  sample  uniformly  each  function  W!^f{x,y)  at  a  rate  proportional  to  2~J . 
Depending  upon  the  choice  of  the  wavelets  4'*(x,}'),  such  a  discretization  can  provide  a  com- 
plete and  orthogonal  representation  [21].  However,  as  in  the  one-dimensional  case,  this  uniform 
sampling  produces  a  representation  which  is  not  translation  invariant.  The  same  translation 
invariance  problem  does  exist  in  non-orthogonal  multiscale  representations  based  on  uniform 
sampling  such  as  the  Laplacian  pyramid  of  Burt  [2]  or  the  DOLP  transform  of  Crowley  [4].  In 
the  next  section,  we  describe  the  two-dimensional  extension  of  the  non-linear  discretization  based 
on  local  maxima.  For  images,  these  local  maxima  provide  the  location  of  the  image  edges. 

5.2.  Wavelet  Transform  Maxima  of  Two-Dimensional  Signals 

We  study  more  particularly  a  two-dimensional  wavelet  decomposition  which  has  only  two 
orientation  selectivities:  a  horizontal  and  a  vertical  one.  This  wavelet  transform  is  defined  with 
two  wavelets  4'Hx,y)  and  4^(j:  ,>')  given  by: 

'i'Kx,y)  =  \\i(x)l,(y)      and       4'2(x.y)  =  ^(x)  \|/(y)  .  (59) 

The  function  \^(x)  satisfies  the  condition  (5)  for  one-dimensional  dyadic  wavelets.  It  is  also 
equal  to  the  derivative  of  a  smoothing  function  Q(x).  The  wavelet  transform  of  a  function 
f{x,y)e  L^(R^)  at  the  scale  2-'  ,  along  the  horizontal  and  vertical  orientations  is  given  by: 

Wif(x,y)  =  /  *  'i'l(x,y)     and    Wlf(x,y)  =  /  *  4'|(A:,y)  .  (60) 

The  condition  (53)  of  two-dimensional  dyadic  wavelet  yields: 

5    |4''(2^co,,2VcOy)l2  +    5    1^2(2^(0;,, 2^ C0y)|2  =  1  .  (61) 

Appendix  6  explains  how  to  choose  the  function  ^(x)  so  that  equation  (61)  holds  and  proves  that 
^(z)  is  a  one-dimensional  smoothing  function.  Let  us  denote  @\x,y)  =  Qix)^(y)  and 
02(j:,)')  =  ^(j:)eCy).  The  functions  Q\x,y)  and  &^(x,y)  are  two  smoothing  functions.  Since 

V|/(x)  =  -^^^^  .the  wavelet  ^\x,y)  and  ^ix,y)  can  be  rewritten: 

nx.y)=  ^Q'j^-y)    and    ^(x,y)=^^^^  .  (62) 

Let  us  denote  by  @lix,y)=  J^  @\^,-^)  and  @lix,y)=  ^  @Hjj-,-^)  ■  Equa- 
tion (60)  and  (62)  yield: 
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Wlf{x,y)  =  2J^(f  *el)ix,y)      and      Wifix,y)  =  2J^(f*@i,){x,y).     (63) 

The  dyadic  wavelet  transform  Wlf(x,y)  and  W^f(x,y)  are  respectively  the  partial  derivative 
along  the  horizontal  and  vertical  directions  of  fix,y)  smoothed  at  the  scale  2-' .  The  local  max- 
ima of  [W^f  (x,y)\  along  x,  for  y  constant,  are  the  points  of  sharper  horizontal  variation  of 
f(x,y)  smoothed  at  the  scale  2-'.  Similarly,  the  local  maxima  of  \Wlf{x,y)\  along  y,  for  x 
constant  are  the  points  of  sharper  vertical  variation.  The  local  maxima  belong  to  curves  in  the 
(x,y)  plane  which  are  the  edges  of  the  image  along  each  direction.  The  detections  of  the  local 
maxima  of  |  Wjjf  (x .y ) |  and  |  W|/ {x,y)\  have  some  similarities  with  a  Canny  edge  detection 
[3].  In  a  Canny  edge  detection  process,  the  smoothing  function  Q(x)  and  ^(x)  are  both  Gaus- 
sians  so  G'CJ^.y)  and  &^{x,y)  are  two-dimensionals  Gaussians.  However,  a  Canny  edge  detccior 
detects  the  local  maxima  of  \Wlf{x,y)\'^+\Wlf(x,y)\'^  at  each  scale  2-'  in  order  to  find  the 
edges  independently  from  their  orientation.  Within  this  wavelet  framework,  we  keep  the  local 
maxima  along  both  the  horizontal  and  vertical  directions  to  discriminate  the  orientation  com- 
ponents of  the  image. 

A  major  problem  in  computer  vision  is  to  imderstand  how  much  information  is  embedded 
in  the  multiscale  edges.  Marr  [24]  made  the  conjecture  that  the  whole  image  can  be  reconstructed 
from  edges.  In  the  following  sections,  we  verify  Marr's  conjecture  by  giving  an  algorithm  that 
reconstructs  the  image  from  the  wavelet  local  maxima  along  the  horizontal  and  vertical  direc- 
tions. 


5  J.  Reconstruction  of  Two-Dimensional  Signals  from  the  Wavelet  Maxima 

The  algorithm  reconstructing  one-dimensional  signals  from  the  maxima  of  their  wavelet 

.,  „ ._    -  ......       '(R")    ^tJ 

^2if(Xyy)y^lfix^y)       7  be  its  dyadic  wavelet  transform.  Lei  us  denote  by  rihesciofall 

sequences  of  functions  g/{x,y) ,  gjHx,y)  ^j.  such  that  for  all  )  e  Z  the  functions  gj\x,y) 
and  Wlf(x,y)  have  the  same  local  maxima  along  x,  at  y  constant,  and  the  functions  g/u,>) 
and  W^f(x,y)  have  the  same  local  maxima  along  y,  at  x  constant.  In  order  to  prove  that  a 
dyadic  wavelet  transform  is  characterized  by  its  local  maxima,  we  must  show  that  there  is  no 
dyadic  wavelet  transform  different  from  the  wavelet  transform  of  f{x,y)  and  which  has  the 
same  maxima.  In  other  words,  the  intersection  of  the  space  V  of  aU  dyadic  wavelet  transform 
with  r  must  be  reduced  to  the  wavelet  transform  of  /  {x  ,y ): 


n 


r  = 


Wif{x,y),Wlfix,y) 


(^^) 


As  in  the  one-dimensional  case,  we  take  a  numerical  approach  to  this  problem  and  develop  an 
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algorithm  for  reconstructing  the  intersection  of  T  and  V.  The  algorithm  also  iterates  on  a  pro- 
jector on  V  and  a  non-linear  projector  on  T.  We  define  a  non-linear  projector  Pr  on  T  which 


transforms  any  sequence  of  functions 


and  that  is  the  closest  to 


gj^ix,y)  ,gj\x,y)     g2  i"^o  ^^  sequence  that  is  in  T 


gjKx,y)  ,gjHx,y)  g2-  ^^  ^^^  '"  section  5.1  that  the  operator 
Pv  =  W  0  W~'  is  an  orthogonal  projection  on  V.  Let  P  =  Py  o  Pp  be  the  composition  of  Pv 
and  Pp.  Any  element  in  the  intersection  of  F  and  V  is  a  fixed  point  of  P.  If  instead  of  defin- 
ing r  only  from  the  maxima  of  the  dyadic  wavelet  transform  of  f  ix,y)  we  use  all  the  extrema. 
one  can  prove  that  F  is  a  convex.  We  can  then  derive  that  the  iteration  on  P  converges  weakly 
to  an  element  in  the  intersection  of  T  and  V  [36].  When  detecting  only  the  maxima  of  the 
wavelet  transform,  we  define  a  non  convex  set  T  so  the  algorithm  is  a  priori  not  guaranteed  to 
converge. 

In  order  to  implement  this  algorithm,  we  first  explain  how  to  define  a  two-dimensional 
dyadic  wavelet  transform  on  a  finite  sequence  of  scales.  The  implementation  of  the  algorithm 
itself  and  the  numerical  results  are  further  described  in  section  6.3. 

6.  Discrete  Model  for  Two-Dimensional  Signals 

6.1.  Discrete  Dyadic  Wavelet  Transform  oflmages 

An  image  is  a  discrete  signal  which  approximates  a  scene  irradiance  at  a  finite  resolution. 
From  such  a  discrete  signal,  we  can  only  define  a  dyadic  wavelet  transfonm  on  a  finite  number  of 
scales.  For  normalization  purposes,  we  assume  that  the  firmest  scale  is  equal  to  1 .  We  suppose 
that  the  wavelet  transform  is  defined  from  the  two  wavelets  given  in  equations  (59).  Let  0(;c  ,_v ) 
be  the  function  whose  Fourier  transform  satisfies: 

|<t>(a),,(0y)|2  =  g  |4''(2^(o,,2^cOy)|2  -H  y  1 4-2(2^(0, ,2V (0,) 1 2  .  (65) 

Let  \\i(x)  be  the  one  dimensional  wavelet  which  defines  ^\x,y)sn(i  ^(x,y)  in  equation  (59), 
Let  ^(x)  be  the  smoothing  function  associated  to  the  wavelet  \^{x)  (equation  (25)).  We  show  in 
appendix  6  that  <^(x  ,y )  can  be  written: 

<D(a:,y)  =  (!)(j: )  ^(y )  .  (66) 

The  function   <I>(x,y)    is  therefore  a  two-dimensional  smoothing  function.  Let  us  denote  by 

<p^(x,y)=  ^X-  '^(-^'-^)-  Let  f(x,y)B  L  (R})  ,  we  denote  by  S2,   tiie  smootiiing  operator 

defined  by: 

S2^fix,y)  =  f  *<^2^ix,y)  . 
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Similarly  to  the  one  dimensional  case,  we  can  prove  that  Sifix,y)   can  be  decomposed  and 
reconstructed  from  the  functions: 


\s2J(x,y),[wlf(x,y)^  j^^^  .  [wlf(x,y) 


lijSJ 


(67) 


This  set  of  functions  is  called  a  finite  scale  two-dimensional  dyadic  wavelet  transform  ot 
Sifix,y). 

At  the  output  of  a  camera  digitizer,  an  image  is  a  finite  energy  two-dimensional  discrete 


signal  D  = 


*n/n 


(n,m)€Z2 


.  With  a  simple  two-dimensional  extension  of  lemma  1  one  can  prove 


that  there  exists  a  function  fix,y)e  L  (R^)  (not  unique)  such  that 

V(«,/n)eZ2    ,    Sif{n,m)  =  dn^   . 


The  discrete  image  can  thus  be  rewritten  D  = 


Sif(n.m) 


(n/n)€Z'' 


(68) 


With  the  class  of  wavelets 


^^(,x,y)  and  ^(x,y)  given  in  appendix  6,  from 


S\f{n,m) 


{n/n)eZ 


J,  we  can  compute  without 


approximation  a  uniform  sampling  of  the  finite-scale  dyadic  wavelet  transform  of  S  \f{x,y) : 


S2J(n,m) 


(«^)6Z''    . 


Wlf(n+w^\m+Wy^) 


(n,m)eZ2 


l<.j<J 


Wlf{n+w^,m+w^') 


(fi  /n  )e  Z'  J   \<j<J 


The  values  {w},w^)  and  {w},w^)  are  constant  shifts  which  define  the  two-dimensional  sam- 
pling grids.  These  shifts  depend  upon  the  choice  of  the  wavelet  4''(^.>')  and  4'2(x,y).  For  any 
scale  2' ,  we  denote 


^h^f  =     V/if{,n+w;t/n+Wy^) 


(n,m)6Z' 


.  W¥f   =     Wlfin+w^^+Wy^) 


,      x.',^  3nd 

(n  ,m  )e  Z- 


sif  = 


S2,f(n,m) 


(n/n)€Z' 


We  call  discrete  dyadic  wavelet  transform  of  the  image  D  the  set  of  two  dimensional  signals 


SU  , 


Wl^f 


\<,j<j  ' 


Wl-'^f 


\<,j<J 


(69) 


For  J  large  enough,  the  samples  of  the  coarse  image  Syf  are  approximately  constant  and  equal 
to  the  mean  of  the  image  D.  Let  us  also  denote  by  VV^  the  operator  which  transforms  a  discrete 
image  D  into  its  discrete  dyadic  wavelet  transform  and  W"''^  the  inverse  operator  which 
reconstructs  an  image  from  a  discrete  dyadic  wavelet  transform.  Appendix  7  describes  two  algo- 
rithms which  implements  both  operators  with  a  complexity  0{N  log{N))  where  A'  is  the 
number  of  pixels  of  the  image  D  .  These  algorithms  are  based  on  separable  convolutions  of  the 
rows  and  colimins  of  the  image  with  one-dimensional  discrete  filters.  The  two  left  columns  of  fig. 
11  give  the  dyadic  wavelet  transform  of  a  grey  circle  in  a  blacic  background.  The  images  along 


the  first  column  correspond  to  the  images  of  the  signals 


W]/f 


\%j<5 


where  as  the  second 


columns  displays 


wm 


\<,j<5 


Page  4 1 
.  The  black,  grey  and  white  pixels  of  these  images  correspond 


respectively  to  negative,  zero  and  positive  sample  values.  These  images  can  be  interpreted 
respectively  as  the  partial  derivatives  along  the  horizontal  and  vertical  directions  of  the  original 
image  smoothed  at  different  scales.  The  image  at  the  top  right  is  S^s/ .  Fig.  12  gives  the  dyadic 
wavelet  transform  of  a  lady  image.  The  coarse  image  S25/  is  approximatively  constant  and  equal 
to  the  average  value  of  the  lady  image. 


Fig.  11-12:  The  top  left  image  is  the  original  image  Sif   and  the  top  right  is  the  smoothed  image 


Sjsf .   The  first  column  from  the  left  gives  the  images     W^i  f 


\<,i<5 


and  the  scale  increases 


from  top  to  bottom.  The  second  columns  is 


w¥f 


J  1^7  25 


.  Black,  grey  and  white  pixels  indi- 


cate respectively  negative,  zero  and  positive  sample  values.  The  third  and  fourth  columns  display 


in  black  the  local  maxima  of  respectively 


Wl^f 


\<,j<S 


and 


Wl^f 


\<j<5' 


In  fig.  11.  the 


darkness  of  maxima  pixels  is  proportional  to  the  amplitude  absolute  value  of  the  corresponding 
maxima.  In  fig.  12,  each  maximum  is  indicated  by  a  black  pixel  independently  from  its  amplitude. 
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Fig.  11:  caption  in  page  41. 
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Fig.  12:  caption  in  page  41. 
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6.2.  Maxima  Detection  of  an  Image  Discrete  Wavelet  Transform 

In  section  5.2,  we  explained  that  we  want  to  detect  at  each  scale  V  the  maxima  ol 
{Wlif  {x,y)\  along  x,  for  y  constant.  These  maxima  belong  to  one-dimensional  curves  in  ihc 
(x,y)  plane.  Because  of  the  discretization,  we  can  only  compute  W^-^/  which  is  a  uniform  sam- 
pling of  W];f{x,y)  at  the  rate  1  along  x  and  y.  Let  us  suppose  for  simplification  that  the  sam- 
pling grids  are  defined  by  w}  =  Wj,'  =  w^  =  w^  =  0.  The  maxima  detection  must  therefore  be  per- 
formed along  each  row   m    of   Wl^'^f    where  we  record  the  maxima  of  the  discrete  one- 


dimensional  signal 


Wlf(n,m)         .   This  detection  is  exact  if  and  only  if  the  maxima  of 

J 

I W^/ (xjn)\  along  x  have  an  integer  abscissa.  Similarly,  the  maxima  of  I  VV|/ {x ,>• ) |  along 
y,  for  X  constant  belong  to  one -dimensional  curves  in  the  (x,y)  plane  but  we  can  only  detect  the 
discrete  maxima  along  each  column    n  of   'W'^'^f .  This  introduces  errors  if  the  maxima  of 

I  ^^/ (''.>')  I    along  y  do  not  have  integer  ordinates.   We  call  maxima  representation  of  the 


VJ]j<^f 


and  along  the  columns  of 


image  D  is  the  set  of  all  maxima  along  the  rows  of 

The  discrete  maxima  detection  is  exact  if  and  only  if  the  maxima  of  |  W^/  {x  jn)\   along  x 
and  the  maxima  of    \y^],f{n,y)\    along  y   appear  at  integer  positions.   Similarly  to  the  one- 


,  plus  the  remaining  coarse  image  S'{,f . 


S\f{n,m) 


{n/n)€7J 


dimensional  case,  we  can  find  a  function  f(x,y)G  L  (R^)  such  that  D  = 

and  whose  finite  scale  wavelet  transform  satisfies  this  property,  only  if  the  function  <P{x  ,y ) 
satisfies  some  further  constraints.  Since  ^(x,y)  =  ^(x)(^(y),  the  results  of  lemma  2  are  easily 
extended.  Appendix  4  proves  that  if  there  exists  a  function  p{x)  which  satisfies  the  conditions 
(34)  of  lemma  2,  then  for  any  discrete  image  D  we  can  find  f{x,y)e  L'CR^)  as  previously 
defined.  One  can  also  show  that  in  this  case,  the  discrete  maxima  along  the  rows  of  W-},-''/  do 
characterize  all  the  maxima  of  Wj,f(_x,y)  along  x,  for  any  ordinate  y  and  the  same  result  is  valid 
for  the  vertical  direction  component.  We  saw  that  the  function  i^{x)  derived  from  the  Haar 
wavelet  admits  a  function  Qix)  satisfying  the  conditions  (34)  whereas  the  cubic  spline  (t)(j: )  ol 
fig.  4(b)  does  noL  This  proves  that  for  wavelets  ^\x,y)  and  ^{x,y)  derived  from  the  Haar 
wavelet,  the  maxima  detection  does  not  introduce  errors  whereas  it  does  if  the  wavelets  wavelets 
^^(x,y)  and  *¥\x,y)  are  obtained  from  the  one-dimensional  cubic  spline  wavelet.  Like  in  one 
dimension,  the  maxima  detection  errors  can  be  viewed  as  the  addition  of  spurious  high  frequen- 
cies in  the  signal  that  we  decompose. 

The  two  columns  on  the  right  of  fig.  1 1  give  the  position  and  amplitude  of  the  local  maxima 
computed  from  the  wavelet  representation  given  by  the  two  columns  on  the  left.  Each  black  or 
grey  point  indicates  the  position  of  a  local  maximum.  The  darkness  of  these  points  is  propor- 
tional to  the  modulus  of  the  maxima  value.  As  expected,  the  local  maxima  of  the  Wj,'^/    and 


Page  45 

^y^f  occur  at  the  border  of  the  circle.  The  amplitude  of  the  maxima  of  'W\;'^f  provide  ihc 
horizontal  component  of  circle  discontinuities  whereas  the  amplitude  of  the  maxima  of  lV|'y 
give  the  vertical  component.  The  two  columns  at  the  right  of  fig.  12  show  the  position  of  ihc 
maxima  of  the  wavelet  transform  given  by  the  left  colunms.  Each  black  pixel  indicates  the  posi- 
tion of  a  local  maximum  independently  of  its  amplitude.  The  diagonal  strucuires  of  the  image 
create  maxima  along  both  the  horizontal  and  vertical  directions.  At  the  finner  scale,  there  are 
many  maxima  for  each  orientation  but  most  of  these  maxima  have  a  negligible  amplitude;  ihey 
are  produced  by  the  image  noise.  The  two  colimins  at  the  left  of  fig.  13  give  for  each  scale  and 
orientation,  the  maxima  of  the  wavelet  transform  of  the  lady  image  whose  amplitudes  are  larger 
than  8.  Many  local  maxima  have  disappeared  with  this  thresholding  and  we  can  see  more  clearly 
the  horizontal  and  vertical  components  of  the  image  edges.  A  proportion  of  68%  of  the  maxima 
are  removed  from  the  maxima  representation  with  this  thresholding  procedure.  In  the  two  right 
columns  of  fig.  13,  the  threshold  is  16  so  we  remove  even  more  maxima  points  from  the 
representation  (79%). 
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Fig.  13:  The  two  left  columns  give  the  position  of  the  maxima  of  the  wavelet  transform  shown  in 
fig.  12  and  whose  amplitude  absolute  value  is  larger  than  8.  Each  black  pixel  indicates  the  posi- 
tion of  a  maximum.  A  proportion  of  68%  of  the  maxima  are  removed  with  this  thresholding.  In 
the  two  right  columns,  we  selected  the  maxima  whose  amplitude  absolute  value  is  larger  than  16. 
The  proportion  of  maxima  removed  is  79%. 
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63.  Numerical  Reconstruction  of  Images  from  the  Wavelet  Maxima 

The  reconstruction  algorithm  is  based  on  the  two  projection  Pv  and  Pp.  The  projection 
Pv  on  the  space  V  of  all  dyadic  wavelet  transform  is  decomposed  into:  Pv  =  W  o  W"'.  The 
discretization  of  the  wavelet  operators  W  and  W~'  was  described  in  the  section  6.1.  For 
discrete  images,  Pv  is  redefined  by: 

P^  =  W  0  W-i^  .  (70) 

If  the  original  image  has  N  pixels,  both  W  and  W"''^  are  implemented  in  0(N  logiN))  so 
the  operator  P^  is  also  implemented  in  0(N  log(N)).  The  other  operator  Pp  is  a  non-linear 
projection  on  the  set  T.  Appendix  5  describes  an  implementation  of  the  corresponding  discrete 
operator  Pj^  with  a  complexity  0(N  log(N)).  The  reconstruction  algorithm  iterates  on  the 
operators  P^  and  Pf^  from  an  initial  sequence  of  discrete  signals  built  with  a  spline  interpola- 
tion between  the  wavelet  maxima. 

The  upper  left  image  of  fig.  14  is  the  original  image  whereas  the  upper  right  image  is  the 
reconstructed  image  from  the  maxima  representation  with  8  iterations.  The  noise  to  signal  ratio 
of  this  reconstruction  is  6.6  10"^.  These  examples  are  computed  with  the  two  wavelets  4^'(i^,y) 
and  ^^x,y)  derived  fi-om  the  one-dimensional  cubic  spline  wavelet.  On  a  video  monitor,  there 
are  no  perceivable  difference  between  the  original  and  the  the  reconstructed  image.  Fig.  15(a) 
gives  the  noise  to  signal  ratio  for  the  reconstruction  of  the  discrete  wavelet  transform  signals 
Wj/f  and     W^'"^/     ,<.<-,,  as  a  function  of  the  number  of  iterations  on  the  operator  V' . 
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For  each  scale  2-' ,  the  noise  to  signal  ratio  does  not  converge  to  zero  but  to  a  minimum  asymp- 
totic error.  The  noise  to  signal  ratio  of  the  remaining  error  decreases  approximatively  by  a  con- 
stant factor  when  the  scale  increases  by  2.  This  proves  that  this  error  is  concentrated  in  the  high 
frequencies.  Like  in  one  dimension,  the  remaining  error  can  be  explained  by  the  spurious  high 
frequencies  introduced  by  the  maxima  detection  procedure.  Indeed,  we  saw  in  the  previous  sec- 
tion that  when  using  the  wavelets  ^\x,y)  and  ^\x,y)  derived  from  the  one-dimensional 
cubic  spline  wavelet,  the  maxima  obtained  from  a  discrete  wavelet  transform  are  not  exactly 
equal  to  the  maxima  of  a  finite  scale  wavelet  transform.  The  corresponding  error  can  be  viewed 
as  spurious  high  frequencies  that  are  added  to  the  signal.  Fig.  15(b)  plots  the  noise  to  signal  ratio 
of  the  image  obtained  by  applying  the  inverse  wavelet  transform  operator  on  the  reconstructed 
discrete  wavelet  transforms.  After  10  iterations  the  noise  to  signal  ratio  is  already  smaller  than 
5.6  10"^.  Because  of  the  remaining  high  frequency  errors,  the  noise  to  signal  ratio  docs  not  con- 
verge to  zero  but  to  a  minimum  value  approximatively  equal  to  2  10"^.  The  remaining  high  fre- 
quency errors  are  well  below  our  visual  sensitivity.  In  practice,  about  10  iterations  are  sufficient 
to  reconstruct  an  image  with  no  perceivable  distortions.  In  order  to  remove  the  high  frequency 
errors  introduced  by  the  maxima  detection,  we  proved  in  the  previous  section  that  we  can  use  the 
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wavelets  ^^{x,y)  and  ^\x,y)  derived  from  the  one-dimensional  Haar  wavelet.  Fig.  16  gives 
the  noise  to  signal  ratio  when  reconstructing  the  lady  image  from  the  wavelet  maxima  obtained 
with  these  wavelets.  After  10  iterations,  the  noise  to  signal  ratio  is  approximatively  3  10"-.  In 
this  case,  the  noise  to  signal  ratio  decreases  steadily  and  reaches  2.5  10"^  after  5000  iterations. 
After  8000  iterations,  the  noise  to  signal  ratio  remains  approximatively  constant  and  equal  to 
10"^.  The  reconstruction  errors  are  concentrated  in  the  high  frequencies  and  the  decay  rate  of  the 
signal  to  noise  ratio  reflects  the  decay  rate  of  the  high  frequency  errors.  As  we  explained  in  the 
one-dimensional  case,  the  reconstruction  of  the  higher  frequencies  is  much  slower  because  the 
reproducing  kernel  operator  P0  is  less  effective  at  the  firmest  scale  2'  of  the  wavelet  transfoirn. 
Our  computations  precision  is  approximatively  10"^.  The  remaining  high  frequency  error  ot 
10"^  is  due  to  a  small  numerical  instability  because  of  the  particular  treatment  of  the  firmer  scale 
by  the  reproducing  kernel  operator.  After  300  iterations,  the  error  made  on  the  value  of  any  pixel 
in  the  reconstructed  image  is  always  smaller  than  0.5.  Since  the  pixel  values  of  the  original 
image  are  coded  with  integers  between  0  and  255,  we  can  recover  exactly  this  image  with  a 
round-off  operation.  The  reconstruction  algorithm  has  been  tested  for  a  large  collection  of  images 
including  special  two-dimensional  functions  such  as  Diracs,  sinusoidal  waves,  step  edges, 
Brownian  noises...  For  all  these  experiments,  the  error  to  signal  ratio  behaves  similariy  to  fig.  15 
or  fig.  16  depending  on  the  wavelets.  Let  us  emphasize  once  more  that  for  image  processing 
applications,  we  need  only  10  iterations  for  reconstructing  an  image  with  no  perceivable  distor- 
tions. 

These  numerical  results  seem  to  indicate  that  the  maxima  representation  is  complete  and 
stable  and  that  the  convergence  errors  are  due  to  the  high  frequencies  errors  introduced  in  the 
maxima  detection  procedure.  We  therefore  conjecture  that  for  a  large  class  of  wavelets  the  max- 


is  a  unique  and  stable  characterization 


ima  of  a  wavelet  transform     W  ^f  ix,y)  ,Wlf  {x  ,y ) 

of  /  (x  ,y ).  This  conjecture  is  based  on  numerical  data  but  is  not  proved  in  this  paper. 

The  stability  of  the  convergence  enables  us  to  slightly  perturbate  the  local  maxima 
representation  and  reconstruct  a  close  image.  The  image  shown  in  the  lower  left  of  fig.  14  is 
reconstructed  from  the  maxima  representation  shown  in  the  two  left  columns  of  fig.  13.  In  this 
case,  we  keep  only  the  maxima  whose  amplitude  absolute  value  is  larger  than  8.  We  thus  remove 
the  maxima  produced  by  the  noise  and  the  light  textures.  As  expected,  the  fine  textures  disap- 
pears in  the  reconstructed  image  but  the  sharp  image  variations  are  not  affected.  The  image 
remains  sharp.  The  image  in  the  lower  right  of  fig.  14  is  a  reconstruction  from  the  maxima 
representation  shown  in  the  right  of  fig.  13.  Even  more  maxima  are  removed  from  the  maxima 
representation  since  the  threshold  on  the  maxima  amplitude  is  twice  larger.  Some  more  textures 
disappear  in  the  reconstructed  image  but  the  strong  edges  and  textures  are  still  unchanged.  This 
type  of  texture  removal  is  similar  to  the  results  obtained  by  Perona  and  Malik  [31]  with 
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anisotropic  diffusion.  For  removing  the  image  noise,  we  can  use  this  reconstruction  technic 
although  we  need  a  more  sophisticated  model  to  select  the  noisy  maxima  that  we  want  to 
suppress.  When  suppressing  local  maxima,  we  also  diminish  the  number  of  values  to  be  coded  in 
the  representation.  This  cleaning  process  thus  also  produces  a  more  compact  representation.  The 
application  to  compact  image  coding  is  further  described  in  section  7.2. 
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Fig.  14:  The  upper  left  image  is  the  original  lady  image.  The  upper  right  image  is  a  reconstruc- 
tion from  the  maxima  representation  shown  in  fig.  12.  This  reconstruction  is  performed  with  8 
iterations  on  the  operator  P^  and  the  noise  to  signal  ratio  is  6.6  10"-^.  The  lower  left  image  is  a 
reconstruction  from  the  maxima  representation  where  the  maxima  amplitude  are  thresholded  by 
8,  as  shown  in  the  two  left  columns  of  fig.  13.  The  lighter  textures  disappear  but  the  image 
remains. sharp.  The  lower  right  image  is  a  reconstruction  with  a  maxima  representation  thres- 
holded by  16,  as  shown  in  the  two  right  columns  of  fig.  13.  Even  more  textures  disappear  while 
the  stronger  edges  and  textures  remain  unchanged. 
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(a)  (b) 

Fig.  15:  (a):  The  ordinate  is  the  noise  to  signal  ratio  when  reconstructing  the  discrete  wavelet 
transform  signals  from  the  wavelet  maxima  of  fig.  12,  on  a  logarithmic  scale.  Each  curve  labeled 
by  a  number  (k)  gives  the  average  of  the  noise  to  signal  ratio  for  the  reconstruction  of  W}/f 
and  W^i'^f .  The  abscissa  is  the  number  of  iterations  on  the  operator  P^ .  These  experiments  are 
done  with  the  wavelets  ^^{x,y)  and  ^\x,y)  associated  to  the  cubic  spline  wavelet  of  fig.  4(a). 
The  noise  to  signal  ratio  does  not  converge  to  zero  because  of  the  high  frequency  errors  intro- 
duced by  the  maxima  detection. 

(b):  The  ordinate  is  the  noise  to  signal  ratio  of  the  image  obtained  from  the  reconstructed  wavelet 
transform.  The  abscissa  is  the  number  of  iterations  on  P*. 
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Fig.  16:  The  ordinate  gives  the  noise  to  signal  ratio  when  reconstructing  the  lady  image  from  the 
wavelet  transform  maxima  obtained  with  the  wavelets  ^^(x,y)  and  ^"^{x.y)  associated  to  the 
one-dimensional  Hoar  wavelet.  The  abscissa  is  the  number  of  iterations  on  the  operator  P^ 
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6.4.  Behavior  of  Maxima  Amplitude  Across  Scales  for  Images 

Equation  (58)  proves  that  the  evolution  of  the  two-dimensional  wavelet  transform  ampli- 
tude characterizes  the  Lipschitz  regularity  of  the  signal  f{x,y).  The  only  difference  with  the 
one-dimensional  situation  is  that  the  wavelet  transform  has  two  orientation  components.  Intui- 
tively the  asymptotic  decay  of     \'W2'f{x,y)\    gives  the  Lipschitz  regularity  of  the  function 
f{,x,y)  for  y  constant  where  as    | W^f  {x,y)\    gives  the  Lipschitz  regularity  of  / (x ,v )   for  x 
constant.   We  can  therefore  extend  the  one-dimensional  procedure  described  in  section  4.4,  in 
order  to  characterize  the  local  shape  of  the  signal  variations  from  the  behavior  of  the  wavelet 
maxima.  Let  us  consider  the  horizontal  component  of  the  wavelet  transform.  We  compute  the 
position  and  amplitude  of  the  extrema  of  W]/  along  each  row.  For  simplicity,  we  can  process 
the  maxima  of  a  row  m  at  different  scales  2-' ,  as  if  it  was  the  maxima  representation  of  a  one- 
dimensional  signal.  We  explained  in  section  4.4  that  from  the  evolution  of  the  maxima  amplitude 
at  the  scales  2',  2^  and  2^,  we  characterize  the  shape  of  the  signal  change  with  three  constants 
K  ,  a  and  Jo-  Let  us  remind  that  K  is  proportional  to  the  amplitude  of  the  signal  change  whereas 
a  and  5o  are  respectively  the  Lipschitz  regularity  (at  the  resolution  1)  and  the  smoothing  scale. 
Equations  (49)  and  (50)  relate  these  constants  to  the  amplitude  of  the  wavelet  maxima  at  the  three 
different  scales.  By  applying  this  processing  to  all  the  maxima  of  W^^f ,  W^/f    and   Wii'^f 
along  the  row  m,  we  characterize  the  variation  of  the  image  intensity  along  the  same  row  m,  in 
the  horizontal  direction.  We  repeat  this  procedure  for  all  rows  m  in  order  to  characterize  the  local 
change  along  the  horizontal  direction  of  the  image  sharper  variation  points.  The  same  computa- 
tions can  be  performed  for  the  maxima  along  the  columns  of  Wl-'^f .  Most  signal  changes  have  a 
horizontal  and  vertical  components  and  thus  create  maxima  along  the  horizontal  and  vertical 
directions.  Let  us  consider  an  edge  which  crosses  a  point  (n  ,m )  of  the  (x,y)  plane  with  an  angle 
G.  This  implies  that  that  the  image  intensity  varies  sharply  when  going  through  the  point  {n  ,m ) 
along  the  direction  of  angle  8  -t-  -y    but  has  no  sharp  variation  along  the  direction  of  angle  9. 
Let  us  suppose  that  this  edge  produces  a  maximum  both  in  Wj'//   and    Wj'"^/   for  1  <  )  <  3,  in 
the  neighborhood  of  {n/n).  To  the  maxima  of  W]j'^f  along  the  row  m  corresponding  to  this 
edge,  we  can  associate  the  constants  K\,  a\  dnA  s\  with  equations  (49)  and  (50)  whereas  we 
obtain  different  constants  Ki,  aa  and  ,^2  from  the  maxima  of  W|-^/   created  by  the  same  edge 
along  the  column  n.  The  Lipschiu  coefficients  ai  and  a:  can  be  interpreted  as  the  Lipschitz 
regularity  of  the  two-dimensional  signal  variation  respectively  along  the  horizontal  and  vertical 
directions.  The  Lipschitz  regularity  a  defined  in  equation  (57)  is  equivalent  to  the  minimum  of 
ai  and  a^.  The  constants  s\  and  S2  give  the  smoothing  scale  along  both  directions.  For  a  step 
edge  of  orientation  9  ;>t  0°  ,  90°.  we  have  ai  =  a2  =  0.  Since  the  wavelet  transform  along  each 
direction  are  partial  derivatives  respectively  along  x  and  y  of  the  image  smoothed,  we  obtain 
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-r7=-  ~  r<3rt(8)  .  For  "comer  points"  where  the  image  varies  sharply  along  any  direction,  the  con- 
stants a\,s\  and  a2  ,  52  can  be  very  different 

As  explained  in  the  one-dimensional  case,  this  simple  analysis  of  the  behavior  of  maxima 
across  scales  is  relevant  only  if  the  maxima  detected  at  the  scale  2^  does  correspond  essentially 
to  one  signal  change  and  not  to  the  aggregation  of  several  small  image  variations.  This  is  a  cm- 
cial  difference  between  the  image  variations  which  are  labeled  "edges"  and  the  one  which  arc 
considered  as  "texture  elements".  This  problem  is  further  discussed  in  section  7.3.  Our  ability  to 
characterize  the  local  properties  of  the  image  edges  opens  a  new  approach  to  edge  detection.  Wc 
argue  in  section  7.3  that  a  fundamental  weakness  of  classical  edge  detectors  is  the  faci  ihai  ihcy 
build  binary  edge  maps  (edge  or  non-edge  point)  which  loose  most  the  image  informaiion  and 
create  instabilities  in  further  processing. 

7.  Applications  of  the  Maxima  Representation  to  Image  Processing 

The  maxima  representation  can  be  viewed  as  a  decomposition  of  the  information  into  mul- 
tiscale  edge  points.  In  the  following  sections,  we  describe  some  applications  of  this  representation 
to  image  coding  and  computer  vision  in  general.  A  key  advantage  of  this  representation  for  classi- 
cal image  processing  is  its  ability  to  decompose  the  image  information  into  meaningful  features. 

7.1.  Maxima  curves 

Along  the  horizontal  direction,  the  local  extrema  of  'W\if{x,y)  along  x,  at  y  constant, 
belong  to  curves  in  the  (x,y)  plane.  However,  since  we  can  estimate  the  position  of  these  extrema 
only  for  y  =m  e  Z,  we  only  detect  isolated  maxima  points.  The  maxima  points  give  a  unifonn 
sampling  along  y  of  the  maxima  curves  of  Wl,f{x,y).  Similariy,  the  maxima  points  of 
W^l/U.y)  that  we  estimate  for  x  =n  €  Z,  give  a  uniform  sampling  along  x  of  the  maxima 
curves  in  the  (x,y)  plane.  For  both  orientations,  the  maxima  curves  are  important  since  ihcy  olicii 
correspond  to  the  boundaries  of  the  image  structures.  In  order  to  rebuild  these  maxima  curve,  wc 
must  regroup  the  maxima  points  that  we  detect.  The  larger  the  scale  1' ,  the  smoother  is  the  van- 
ation  of  the  maxima  amplitude  along  each  maxima  curve.  Indeed,  we  saw  that  Wlfix,y)  and 
Wlf(x,y)  are  the  partial  derivatives  of  f(x,y)  smoothed  at  the  scale  V .  We  are  therefore 
going  to  group  together  maxima  points  which  are  neighbors  and  have  a  close  amplitude.  A  local 
maximum  detected  along  a  row  r  of  Wjj'^f  is  grouped  with  a  maximum  along  a  row  r-i- 1  if  both 
their  position  and  amplitude  are  close.  Each  maximum  is  grouped  at  most  with  one  maximum  in 
the  previous  row  and  one  maximum  in  the  next  row.  Similariy,  we  group  together  the  maximum 
points  of  Wl'^f  detected  along  consecutive  columns. 
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The  length  of  a  maxima  curve  and  the  average  of  the  maxima  amplitude  along  the  curve  arc 
two  important  descriptors.  Let  us  call  maxima  curve  amplitude  the  multiplication  of  these  two 
values  which  is  equal  to  the  sum  of  the  maxima  amplitude  along  the  maxima  curve.  The  maxima 
curves  whose  amplitudes  are  large  in  absolute  value  are  the  long  curves  or  the  one  having  high 
amplitude  maxima.  Instead  of  thresholding  the  maxima  representation  based  on  the  amplitude  of 
individual  maxima  as  it  was  done  in  fig.  13,  one  can  make  a  thresholding  of  the  maxima  curves 
based  on  their  amplitude.  The  image  at  the  top  left  of  fig.  17  is  the  original  image  whereas  at  the 
top  right  is  a  reconstruction  from  a  maxima  representation  where  we  only  kept  the  maxima  curves 
whose  amplitude  are  larger  than  16.  The  four  images  at  the  bonom  left  display  all  the  maxima  ol 
the  wavelet  representation  at  the  scales  2'  and  2^,  along  the  horizontal  and  vertical  orientations. 
The  images  at  the  bottom  right  show  the  maxima  curves  whose  amplitude  absolute  value  arc 
larger  than  16,  for  the  same  orientations  and  scales.  This  simple  thresholding  procedure  selects 
most  of  the  image  important  edges.  Although  69%  of  the  maxima  points  have  been  removed 
from  the  original  representation,  few  differences  do  appear  between  the  original  image  and  the 
reconstructed  one.  In  the  following  we  consider  the  maxima  curves  as  the  elementary  features  ol 
the  maxima  representation. 
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Fig.  17:  The  original  image  is  at  the  top  left.  The  top  right  image  is  reconstructed  from  a  max- 
ima representation  where  at  all  scales,  we  only  keep  the  maxima  curves  whose  amplitude  are 
larger  than  16.  The  four  images  at  the  bottom  left  display  all  the  maxima  of  the  wavelet 
representation  at  the  scales  2'  and  2^,  along  the  horizontal  and  vertical  orientations.  The 
images  at  the  bottom  right  show  the  maxima  curves  whose  amplitude  absolute  value  are  larger 
than  16,  for  the  same  orientations  and  scales. 
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7.2.  Compact  Image  Coding 

Image  coding  is  a  fundamental  problem  not  only  in  classical  image  processing  but  also  in 
computer  vision.  Any  high  level  description  of  the  image  can  also  be  viewed  as  a  coding.  An  effi- 
cient procedure  to  build  a  compact  image  coding  is  to  group  similar  features.  It  is  generally  more 
economical  to  code  similar  features  as  a  group  rather  than  each  feature  independently.  Within 
each  group,  the  features  might  not  be  identical  which  requires  to  record  their  differences  like  in 
classical  predictive  coding  technics.  The  principles  of  grouping  in  perception  have  been  studied 
and  formalized  by  the  Gestalt  psychology  in  the  1930's.  In  computer  vision,  Marr  [24]  also  pro- 
posed to  use  percepnjal  grouping  to  build  a  more  effective  image  representation  from  his  "prima! 
sketch".  More  recently  Lowe  [19]  made  a  thorough  analysis  of  the  grouping  approach  to  image 
analysis.  In  section  7.1,  we  implemented  a  first  grouping  step  by  chaining  together  the  individual 
local  maxima  to  recover  the  local  maxima  curves.  More  sophisticated  high  level  grouping  is  also 
possible. 

One  approach  to  compact  image  coding  is  to  divide  into  three  groups  the  image  information 
provided  by  the  maxima  curves.  The  first  group  is  the  image  "noise"  which  is  defined  as  the 
information  that  can  be  removed  without  affecting  the  image  "quality".  Of  course  the  quality  cri- 
teria depends  upon  the  type  of  information  that  we  want  to  obtain  from  the  image.  With  a  simple 
thresholding  operation  on  the  maxima  curve  amplitude  we  reconstructed  in  fig.  17  an  image  that 
is  very  close  from  the  original  one  from  a  perception  point  of  view  although  we  removed  69%  oT 
the  maxima  curves  from  the  image.  This  diminishes  considerably  the  number  of  bits  for  coding 
the  image.  A  more  thorough  study  is  needed  in  order  to  understand  what  maxima  curves  can  be 
removed  without  affecting  the  visualization  of  the  reconstructed  image.  The  second  categon,' 
regroups  the  "textures"  which  are  defined  for  coding  purposes  as  the  image  variations  that  do  not 
have  an  information  content  by  themselves  but  as  a  group.  For  example,  we  do  not  consider  each 
hair  of  the  lady  shown  in  fig.  17  but  the  hairs  as  a  whole.  Texture  discrimination  is  one  of  the 
most  difficult  problem  that  is  further  discussed  in  the  next  paragraph.  One  can  observe  thai 
important  distortions  can  be  introduced  in  textures  without  affection  tiieir  visualization.  Textures 
can  thus  be  coded  with  relatively  few  bits.  Finally  the  third  group  are  the  "edges"  which 
correspond  to  the  maxima  curves  that  must  be  coded  individually  because  they  are  isolated  and 
are  significant  by  themselves.  Edges  cannot  be  modified  substantially  without  producing  per- 
ceivable changes  in  the  image.  The  edge  maxima  curves  can  be  coded  efficienUy  with  B-splincs. 

In  order  to  make  such  an  image  coding,  we  must  discriminate  the  maxima  curves  that  arc 
considered  as  noise,  texture  or  edge.  One  can  use  the  geometrical  properties  of  the  maxima 
curves  for  this  purpose.  The  noise  generally  corresponds  to  the  maxima  curves  that  are  short  and 
not  regular.  On  the  contrary,  long  and  smooth  maxima  curves  are  often  important  boundaries  of 
some  image  structure.  The  behavior  across  scales  of  the  maxima  curves  is  another  important  cue 
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to  discriminate  textures  from  edges.  We  mentioned  in  section  7.4  that  textures  generally  create 
maxima  at  a  fine  scale  2-'  that  breaks  into  several  others  maxima  when  the  scale  decreases.  \Vc 
do  believe  that  the  wavelet  maxima  representation  is  particularly  well  adapted  for  high  compres- 
sion image  coding  because  of  its  ability  to  discriminate  the  different  type  of  features  that  arc 
important  from  a  perception  point  of  view. 

7.3.  Edge  Detection,  Texture  discrimination  and  Regularization 

There  is  no  need  to  further  motivate  the  importance  of  edges  for  image  pattern  recognition 
since  many  existing  vision  algorithms  are  based  on  sharp  variation  detection.  A  fundamental 
application  of  edge  detection  is  to  find  the  boundaries  of  the  image  structures  and  segment  the 
image  into  "meaningful"  regions.  After  many  years  of  research  on  trying  to  find  the  "optimal" 
local  edge  detector,  more  recentiy  people  have  tried  to  include  in  the  edge  detection  process  some 
a  priori  knowledge  on  the  image.  This  leads  to  a  regularization  formulation  of  edge  detection.  A 
regularization  method  consists  of  using  some  a  priori  knowledge  on  the  expected  solution  of  an 
ill-defined  problem  in  order  to  compute  an  admissible  solution.  Most  problems  in  low-level 
vision  are  ill-defined  and  need  to  be  stabilized  with  a  regularization  procedure  [32].  The  classical 
regularization  procedure  make  a  smoothness  assumption  on  the  image  intensity  which  is  not  valid 
in  presence  of  textures.  We  discuss  the  wavelet  transform  approach  for  regularizing  non-smooih 
problem  by  using  assumptions  on  image  singularities. 

To  regularize  edge  detection,  Mumford  &  Shah  [30]  as  well  as  Blake  &  Zis.scnman  fl] 
associate  edge  finding  to  the  solution  of  a  variational  problem.  The  image  is  considered  as  a 
piecewise  smooth  surface  with  noise.  The  regularization  process  tries  to  recover  the  "noisc-frcc" 
image  by  minimizing  an  energy  functional  which  expresses  the  a  priori  knowledge  on  the  image 
smoothness  as  well  as  a  constraint  on  the  total  length  of  the  curves  along  which  the  image  is 
singular.  These  curves  are  the  image  edges.  This  energy  can  also  be  interpreted  as  the  Gibbs 
energy  of  a  Markov  random  field  [9].  The  minimum  of  the  energy  functional  can  thus  also  be 
interpreted  as  the  maximum  of  the  a  posteriori  probability  distribution  of  the  Markov  field,  given 
the  noisy  data  provided  by  the  image.  The  minimization  of  this  energy  functional  is  a  difficult 
problem  [29]  and  requires  expensive  numerical  algorithms  such  as  simulated  annealing  [9].  The 
major  inconvenience  of  this  regularization  scheme  is  that  the  global  smoothness  assumption  is 
not  valid  for  most  images.  Indeed,  images  generally  include  irreg<jlar  textures  and  are  thus  non- 
smooth  almost  everywhere.  The  Mumford  &  Shah  regularization  model  of  edge  detection  has  its 
roots  in  the  belief  that  an  edge  detection  process  is  ultimately  a  binary  process  which  classifies 
the  image  points  into  edges  and  non-edges.  The  underlying  assumption  is  that  there  arc  few 
singularities  in  the  image  and  that  all  these  singularities  are  discontinuities  (step  edges).  For 
example,  this  model  does  not  distinguish  an  occlusion  edge  from  a  much  smoother  edge  o(  :i 
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shadow.  It  will  also  classify  all  the  irregularities  of  textures  either  into  discontinuity  points  or 
smooth  points.  This  gives  a  very  crude  description  of  the  image  information  which  often  leads  to 
unstable  recognition  process  at  latter  stages.  One  might  argue  that  an  important  advantage  of  a 
process  which  only  labels  points  as  edge  or  non-edge  is  to  considerably  reduce  the  amount  of  data 
representing  the  image  and  thus  speed  up  the  following  computations.  We  think  that  in  most 
situations,  this  view  is  inexact.  In  order  to  compensate  the  lack  of  low-level  infonnation  provided 
by  a  classical  edge  detector,  most  vision  systems  have  to  use  higher  level  cues  which  are  often 
ad-hoc  and  requires  much  more  heavy  processings. 

We  believe  that  edge  detection  should  not  only  separate  smooth  regions  from  discontinui- 
ties but  rather  characterize  the  type  of  discontinuities  which  appear  in  the  image.  Even  if  ihc 
image  intensity  is  extremely  irregular  almost  everywhere,  there  are  some  singularities  thai  one 
might  want  to  distinguish  from  others  because  of  their  local  shape  or  their  geometrical  propcnics 
in  the  (x,y)  plane.  For  example,  in  the  lady  image  shown  in  fig.  12  the  irregularities  on  the  fur 
texture  on  the  hat  are  very  different  from  the  intensity  discontinuities  at  the  border  of  the  lady's 
shoulder.  The  fur  is  an  important  aspect  of  the  image  and  cannot  be  considered  as  noise.  The  rcg- 
ularization  process  should  then  try  to  differentiate  these  different  singularities  based  on  their  type 
and  geometrical  properties.  We  saw  in  section  6.4  that  one  can  characterize  the  edge  types  from 
the  decay  of  the  wavelet  maxima  across  scales.  We  can  recover  the  Lipschitz  regularity  of  these 
sharp  variation  points  and  detect  any  smoothing  in  these  points.  This  regularization  approach  can 
be  applied  for  regularizing  non-smooth  images.  Instead  of  making  a  global  smoothness  assump- 
tion, one  can  impose  some  constraint  on  the  type  of  singularities  appearing  in  the  solution  of  an 
ill-defined  problem. 

One  of  the  most  challenging  problems  in  low-level  vision  is  texture  discrimination.  A  tex- 
ture is  a  region  where  the  image  intensity  might  vary  considerably  but  which  is  perceived  as  uni- 
form because  the  variations  look  alike.  Right  now,  the  human  visual  system  provides  the  only 
reliable  measurements  of  texture  homogeneity.  Texture  discrimination  is  thus  still  an  experimen- 
tal problem  where  the  only  measurement  of  success  is  the  comparison  against  the  human  visual 
system.  Texture  discrimination  and  edge  detection  are  two  aspects  of  the  same  problem.  .'\n 
image  sharp  variation  might  be  considered  as  a  texture  element  at  one  scale  or  as  an  important 
boundary  (edge)  at  a  firmer  scale.  The  scale  is  thus  a  particulariy  important  parameter  for  texture 
analysis.  When  discriminating  textures,  one  looks  for  undedying  regularities  in  an  irregular 
region.  The  ability  of  the  maxima  representation  to  characterize  the  irregularity  types  is  important 
for  this  problem.  One  of  the  difficulties  with  textures  that  we  mentioned  in  section  4.4  is  that  the 
image  variations  are  not  isolated.  We  thus  need  to  develop  a  firmer  model  to  interpret  the  interac- 
tion of  the  signal  variations  on  the  maxima  amplitude  at  the  larger  scales.  The  reconstruction 
algorithm  from  maxima  provides  an  interesting  tool  for  experimental  studies  on  textures.  .Any 
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model  can  be  tested  by  synthesizing  textures  with  the  reconstruction  algorithm  and  by  comparing 
the  results  of  the  model  against  the  performance  of  a  human  observer. 

7.4.  Pattern  Matching 

Multiscale  representations  have  already  found  many  applications  for  matching  signals. 
Indeed,  pattern  matching  algorithms  can  be  implemented  efficiently  with  coarse  to  fine  pro- 
cedures. At  a  large  scale,  there  are  few  edges  in  the  images  and  the  edges  describe  the  larger 
strucuires.  These  structures  are  generally  less  sensitive  to  the  distortions  which  might  exist 
between  the  patterns  that  should  be  matched.  The  matching  is  quickly  performed  at  large  scales 
because  of  the  small  number  of  edges.  We  then  refine  progressively  the  scale  and  use  the  infor- 
mation provided  by  matches  at  larger  scales  to  speed  up  the  matching  process.  Stereo  matching 
algorithms  based  on  this  principle  have  been  developed  by  Marr  &  Poggio  [26],  Crimson  [I0| 
and  others.  For  matching  signals,  we  compare  the  multiscale  edges  which  requires  an  undcriying 
distance.  In  most  existing  algorithms,  the  distance  is  ad-hoc.  The  wavelet  model  enables  us  to 
build  a  distance  on  multiscale  edges  by  using  the  amplitude  of  the  wavelet  maxima.  By  interpo- 
lating the  wavelet  maxima  with  dilated  spline,  we  can  define  a  simple  L  (R^)  metric  which 
approximates  the  mean  square  distance  between  the  original  wavelet  transforms.  An  algorithm 
for  matching  one-dimensional  signals  has  already  been  implemented  from  zero-crossing  (instead 
of  maxima)  with  a  similar  idea  [20].  The  completeness  of  the  wavelet  multiscale  edges  show  thai 
such  a  matching  procedure  does  not  need  to  be  refined  by  a  local  correlation  between  pattern  pix- 
els as  it  is  generally  done  in  most  algorithms.  In  stereo  matching,  the  patterns  are  translated 
along  the  horizontal  direction  (after  rectification)  but  not  rotated.  It  is  thus  natural  to  decompose 
the  image  information  along  the  horizontal  and  vertical  directions.  However,  for  some  other  type 
of  matching  like  motion  correspondence,  the  matching  must  also  recover  rotation  parameters.  Wc 
can  then  either  decompose  the  signal  along  more  orientations  in  order  to  be  almost  rotational ly 
invariant  or  combine  together  the  horizontal  and  vertical  orientations.  Intuitively,  since  ihc 
wavelet  transform  computes  the  gradient  along  the  horizontal  and  vertical  directions,  wc  can 
recover  the  gradient  along  any  direction  (at  each  scale  2-').  A  simple  model  for  this  integration 
was  mentioned  in  section  6.4. 

8.  Conclusion 

We  saw  that  multiscale  edges  can  be  obtained  from  the  local  maxima  of  a  wavelet 
transform  and  we  explained  how  to  reconstruct  the  original  signal  from  these  local  maxima.  The 
reconstruction  algorithm  can  restore  exactly  the  image  and  thus  gives  a  verification  of  Marr  con- 
jecture on  the  completeness  of  multiscale  edges.  However,  the  mathematical  study  of  the  com- 
pleteness problem  remains  on  open  question.  We  also  proved  that  one  can  characterize  the  type 
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of  the  edge  points  from  the  decay  of  the  wavelet  maxima  amplitude  across  scales.  This  is  particu- 
larly important  to  discriminate  the  different  image  sharp  variations.  The  wavelet  ma.xima 
representation  is  a  new  reorganization  of  the  image  information  in  order  to  develop  solutions  to 
image  processing  problems,  from  the  properties  of  the  image  edges.  Although  we  emphasized 
image  processing,  this  representation  has  potential  applications  for  processing  other  types  of  sig- 
nals. Its  ability  to  characterize  the  signal  irregularities  and  its  translation  invariance  arc  tv\o 
important  properties  for  the  analysis  of  transient  signals  in  general. 
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Appendix  1 

Proof  of  Lemma  1 


For  any  finite  energy  discrete  signal  D  = 


neZ 


,  we  want  to  find  /(x)  e  L^(R)  such 


that 


V/i  €  Z    .    5i/(«)  =  dn  .  (73) 

Let  f(x)€  L  (R),  by  definition  we  have  S\f(n)=f  *  (t>(/i)    .  This  convolution  product  can  be 
rewritten  as  an  inner  product  in  L  (R) :  5i/(n)  =  </U) ,  <\)(n-x)>  .  Let  U  be  the  vector  space 


generated  by  the  family  of  functions 


i^ix-n) 


neZ  ■ 


If  this  family  is  a  basis  of  U  then  for  any 


discrete  sequence 


neZ 


of  finite  energy,  there  exists  f{x)€  L'CR)  satisfying  equation 


(73).  One  can  show    [22]  that  the  family 
Fourier  transform  $((o)  satisfies: 


<^ix-n) 


nez   is  a  Hilbert  basis  if  and  only  if  the 


The  values 


3Ci  >0  ,  3C2>0  ,  Vcoe  R  ,     Ci  <     ^    |$(0>f2/zJt)|2  <  C2 

2, 


^^2  characterize  the  orthogonal  projection  of  fix)e  L  (R)   on   U  .  This 

orthogonal  projection  can  be  interpreted  as  an  approximation  at  the  resolution   1  of  the  function 
f(x)  [21]. 
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A  Particular  Class  of  One-Dimensional  Dyadic  Wavelets 

This  appendix  defines  the  class  of  wavelets  used  for  implementation  of  discrete  algorithms. 


From 


Sxf(n) 


neZ 


we  want  to  be  able  to  compute 


Syfin) 


neZ  ' 


W2,fin+w) 


neZ 


\<J<J 


with  discrete  convolutions,  for  some  constant  w.   If  7  =  1,  this  implies  that  we  can  compute 


S2./(n) 


neZ 


by  convolving 


Sxf{n) 


n€  Z 


with  a  discrete  filter    H.  In  other  words  the 


Fourier  series  of 


Stfin) 


neZ 


is  equal  to  the  Fourier  series  of 


5i/(n) 


/1€Z 


multiplied  by  a 


2jt  periodic  function  //((o) .  The  Fourier  series  of  these  two  signals  are  respectively 

5  /  *<t>(i )«""""    and      2  /  *'!>2(«)e"'"'^  • 
By  applying  the  Poisson  formula,  we  can  rewrite  these  two  series: 


(74) 
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X  / (CW-2/I  Ji)  $(o>+2«  n)    and        X  /(o>+2/in)<^(2(a>i-2n7i))    .  (75) 


The  left  series  is  equal  to  the  right  series  multiplied  by  //((o)  for  all  /  (co),  if  and  only  if 

(f(2co)  =  //((o)(^(co)  .  (76) 

Since  |  (f  (0)  |  =  1 ,  we  must  have  |  //  (0)  |  =1  .  If  we  cascade  equation  (76),  we  obtain  a  necessary 
condition  on  ^((o) : 

$((0)  =  q//(2-P(0)  .  (77) 

Conversely,  if  the  2k  periodic  function  //(co)  satisfies 

|//((0)|2  +    |//((i>frt)|2  <   1    .  -         (78) 

then  one  can  show  [22]  that  the  function  (\)ix)  whose  Fourier  transform  is  defined  by  equation 
(77)  is  a  function  in  L^(R).  The  function  //(co)  can  be  interpreted  as  the  transfer  function  of  a 
discrete  low-pass  filter. 

Let  us  now  characterize  the  corresponding  wavelet  \\tix)  .  As  a  consequence  of  equation 
(29),  we  have 

IV(2co)|2  =  |$({o)|2-  |(f(2a))|2  .                                       (79) 

Substituting  equation  (76)  in  equation  (79)  yields 

\i}(2co)  =  G  ((0)  $(co)  .with  (80) 

|G((D)|2  +  |//(co)|2  =  1  .  (81) 

The  function  G  (co)  is  chosen  2k  periodic  and  can  be  interpreted  as  the  transfer  function  of  a 
high-pass  filter.  If  the  inequality  (78)  is  an  equality  then  we  can  choose  G((o)  =  e~"^H {(i>+k)  . 
One  can  show  that  this  defines  a  wavelet  V|/(j:)  which  generates  a  wavelet  onhonormal  basis  as 
described  in  section  2.4  [21].  However,  this  condition  is  restrictive  and  does  not  bring  any  advan- 
tage in  the  local  maxima  context. 

For  the  maxima  model,  we  want  to  build  a  wavelet  \\iix)  equal  to  a  first  order  derivative  of 
a  smoothing  function  6(x).  This  implies  that  y((o)  must  have  a  zero  of  order  1  in  (0  =  0.  Since 
1 0(0)  1=1,  equation  (80)  yields  that  G  (co)  must  only  have  a  zero  of  order  1  in  co  =  0.  The  func- 
tions G(co)  and  //(co)  are  related  by  equation  (81),  this  condition  is  equivalent  to 

//(CO)  =  1  +  0(co2)    ,  (82) 

for  0)  in  the  neighborhood  of  zero.  For  most  applications  in  image  processing  it  is  important  to 
use  a  wavelet  \^ix)  which  is  antisymmetrical  in  order  to  operate  similarly  on  the  left  and  the 
right  neighborhood  of  a  point.  This  constraint  is  satisfied  if  the  function   //(w)    is  real  and 
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symmetrical  around  zero  and  G(co)  is  imaginary  pure  and  antisymmetrical  around  zero.  Given 
all  these  constraints,  there  still  is  a  large  possible  choice  for  //(co).  Daubechies  [7]  showed  how 
the  choice  of  functions  //(co)  and  G(co)  influences  the  smoothness  and  decay  rate  at  infinity  of 


\|/(x).    Let   H  = 


ngZ 


and    G  = 


gn 


neZ 


denote  respectively  the  discrete  filters  whose 


transfer  function  are  respectively  //((O)  and  G((0).  The  function  //(co)  corresponding  to  the 
function  ^(x)  and  y(j:)  shown  in  fig.  4  is 


//(CO)  =  (coj(-^))4  . 


(83) 


The  impulse  response  of  the  corresponding  filters  H  and  G  are  given  in  table  2.  The  filter  H  is 
symmetrical  around  zero  and  has  only  three  non-zero  coefficients.  Tlie  filter  G  is  antisymmetri- 
cal around  zero  and  its  coefficients  are  exponentially  decreasing.  We  only  give  the  first  5  coeffi- 
cients. The  corresponding  wavelet  \\i{x)  is  a  cubic  spline  having  an  exponential  decay. 


n 

h„ 

gn 

0 

0.3750 

0.00000 

1 

0.2500 

0.59261 

2 

0.0625 

0.10872 

3 

0.0000 

0.01643 

4 

0.0000 

0.00008 

Table  2:  First  five  coefficients  of  the  impulse  response  of  the  filters  H   and  G  corresponding  to 
the  cubic  spline  wavelet  in  fig.  4(a). 


The  function  //(co)  corresponding  to  the  Haar  wavelet  (equation  (36))  is  given  by 

//(CO)  =  e    ■^cos(-|-)  . 

The  impulse  responses  of  the  corresponding  filter  H  and  G  are 

_\  \  if  n  =0  OT  n  =  1 
"n  ~  1  0  otherwise 


(84) 


gn    =1 


1    if   Al  =0 

-1  if  n  =  1 
0  otherwise 


(85) 


(86) 


For  computing  the  Lipschitz  regularity  a  and  the  smoothing  scale  sq  in  the  numerical 
experiments  of  section  4.4,  we  used  a  wavelet  y(x)  characterized  by: 
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H((ii)  =  e  '^  (cos(^))'^  . 

The  first  3  coefficients  of  the  impulse  response  of  the  filters  H  and  C  are  given  in  tabic  3.  The 
coefficients  of  these  filters  are  respectively  symmetrical  and  antisymmetrical  around  -^.  The 
corresponding  wavelet  is  a  quadratic  spline.  It  is  continuously  differentiable  and  antisymmeirical 
around    4-.  The  primitive    Qix)   of  this  wavelet  satisfies  approximatively  equation  (38)    for 


p(2^5o)  =  ^/22J+io^ 


n 

K 

gn 

1 

2 
3 

0.3750 
0.1250 
0.0000 

0.5798 
0.0869 
0.0061 

Tah\e  3:  First  three  coefficients  of  the  impulse  response  of  the  filters  H   and  G.  The  impulse 
response  of  these  filters  is  respectively  symmetrical  and  antisymmetrical  around  y . 


Appendix  3 

Fast  Wavelet  Algorithms  for  One-Dimensional  Signals 

This  appendix  describes  an  algorithm  for  computing  a  discrete  wavelet  transform  and  the 
inverse  algorithm  that  reconstructs  the  original  signal  from  its  wavelet  transform.  We  suppose 
that  the  wavelet  v(x)  is  characterized  by  the  two  discrete  filters  H  and  G  described  in  appen- 
dix 2.  We  denote  Hp  and  Gp  the  discrete  filters  obtained  by  putting  2^  -  1  zeros  between  each 
coefficients  of  the  filters  H  and  G.  The  transfer  function  of  these  filters  is  respectively  //(2''co) 
and  Gi2P(o).  We  also  denote  by  Hp  and  Gp  the  filters  whose  transfer  function  are  respec- 
tively 7/(2'' (o)  and  G  (2^  CO)  (complex  conjugates  of //(2''a))  and  G(2Pco)).  We  denote  by 
A  *B   the  convolution  of  two  discrete  signals  A  and  B. 

The  following  algorithm  computes  the  discrete  wavelet  transform  of  the  discrete  signal 
Sif .  At  each  scale  2J ,  it  decomposes  5^,^  into  S^f  and  W^f . 
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j  =  0 

while  (j  <  J) 

Wi4  =  Sif  *  Gj 

Si4  =  Sif  *  Hj 

end  of  while 

The  proof  of  this  algorithm  is  based  on  the  properties  of  the  wavelet  \(/(a:)  described  in 


appendix  2.   If  the  original  signal 


5i/(/i) 


ns  Z 


has  A'   non-zero  samples,  then  each  signal 


S^f  and  W^f  has  A^  non-zero  sample.  Since  there  are  at  most  log{N)  scales,  the  complexity 
of  the  algorithm  is  OCA'  log{N))  .  The  constant  depends  upon  the  number  of  non-zero  coeffi- 
cients in  the  filters  H  and  G .  In  order  to  avoid  border  problems,  one  can  extend  the  signal  by 
symmetry  around  the  first  and  the  last  sample.  If  H  is  symmetrical  around  0  and  G  antisym- 
metrical  around  0,  each  signal  S^f  and  W^f  are  resp)ectively  symmetrical  and  antisymmetri- 
cal  around  the  first  and  last  sample.  For  the  filter  H  and  C  of  the  Haar  wavelet  the  symmetry 
is  slightly  different  due  to  a  shift  of  -L-  in  the  symmetry.  Indeed,  the  shifting  constant  w  that 
defines  the  sampling  grid  of  W^f  depends  upon  the  phase  of  the  transfer  function  G  (co)  of  the 


filter  G.  For  a  Haar  wavelet,  w  = 


_  1 


T- 


wif  = 


W^fin+].) 


ne  Z 


The  inverse  wavelet  transform  algorithm  reconstructs  Sif  from  the  discrete  dyadic 
wavelet  transform.  At  each  scale  2-' ,  it  reconstructs  S^.J  from  5^/  and  W^f  .  The  com- 
plexity of  this  reconstruction  algorithm  is  also  OiN  logiN)). 


while  (J  >  0) 

Si4  =  Wif  *  G,-x  +  Sif  *  fij-i 

end  of  for 
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Appendix  4 

Proof  of  Lemma  2  and  Two-Dimensional  Extension 

In  this  appendix  we  study  the  errors  due  to  the  discretization  of  the  finite  scale  wavelet 
transform  when  estimating  the  position  and  amplitude  of  the  maxima.  We  prove  lemma  2  and 
make  a  simple  two-dimensional  extension  of  these  results. 

At  each  scale  2-',  l<j<J,we  must  compute  the  maxima  of  W2jf{x)  from  its  uniform 
sampling  W^f .  For  this  purpose,  we  must  insure  that  the  maxima  of  Wj/ZU)  have  integer 
abscissa.  Let  us  suppose  that  there  exists  5ip(a:)  =  p  *  (t>(j:)  such  that 


SipCx)  =  S]p{-x)     Vx  €  R 

5ip(x)  =  0     if    Ul  >1 

— J         <  0    for   0  <  j:  <  1 
ax 

Sipix)  +   Sip{\-x)  =1     if  0<x  <  1 


(87) 


The  function  S  ipU)  is  a  non-oscillatory  interpolation  of  a  discrete  Dirac 

J  1  if  n  =0 
S\Pin)  =  |o  if  Ai  *0  • 

By  non-oscillatory  we  mean  that  it  has  no  maxima  between  the  samples  with  integer  abscissa. 


Let  D  = 


dn 


1€Z 


and  let  us  define 


fix)  =    ^  d„  p{x-n) 


(88) 


This  equation  yields  that 


S^f{x)  =    X  d„Sipix-n) 


(89) 


We  can  easily  prove  that  S\f(x)  is  a  non-oscillatory  interpolation  of  the  samples  of  D.  Let  us 
show  that  for  the  class  of  wavelets  described  in  appendix  2,  for  all  scales  2^  >  1 ,  we  also  have 

W2^f(x)  =    f^W-^fin)  Sip{x-n)  .  (90) 

Equations  (77)  and  (80)  yield  that  the  Fourier  transform  of  a  wavelet  \^^ix)  satisfies  for;  >  0 

^{2J(ii)  =  M;(co)$(a))  , 

where  Mj(co)  =  G(2>-i(o)  "n  ^(2"''w)  is  a  In  periodic  function.  Let  Mj   be  the  discrete  filter 
whose  Foiuier  series  is  Mj  (co).  We  have 

W2,f(x)=f  *  x|/2;(x)=/  *  (Mj  *  (t))(x)   so 
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W2,f(x)  =  Mj  *SJix)  . 
From  this  equation  and  equation  (89),  we  derive  that  H^2^/(x)  can  be  written 

W^f(x)  =  J^w„  Sip(x-n)  . 

Since  S\p(x)  satisfies  the  conditions  (87),  we  obtain  w„  =W2,f{n).  Equation  (90)  yields  that 
the  maxima  of  W-^fix)  have  integer  abscissa  and  are  therefore  identical  to  the  discrete  maxima 
ofW^f.  This  proves  that  if  there  exists  p(x)  which  satisfies  the  conditions  (34)  of  lemma  2 
(conditions  (87))  then  for  any  discrete  signal  D,  the  function  Sif(,x)  given  in  (89)  interpolates 
the  samples  of  D  and  the  maxima  of  Wjifix)  have  integer  abscissa.  This  concludes  the  proof 
of  lemma  2. 

The  maxima  representation  is  defined  from  the  amplitude  and  position  of  the  discrete  max- 
ima of  W^f  ,J  <j<\.  We  suppose  implicitly  that  the  discrete  maxima  of  W^f  are  identical 
to  the  maxima  o{W2,fix).  For  any  dyadic  wavelet  \\iix),  it  is  not  true  in  general  that  there  exists 
a  function  f{x)e  L^(R)  such  that  SJix)  interpolates  the  samples  of  D  and  the  maxima  of 
^2if(.^)  have  integer  abscissa.  Let  us  suppose  for  example  that  the  function  (^(x)  associated  to 
\\i(x)  has  a  support  larger  than  2.  Lemma  1  proves  that  there  exists  p(x)  e  L^(R)  such  that 
S  \p(x )  interpolates  a  discrete  Dirac: 

„     ,  ,      J  1  if  /I  =0 
S\Pin)  =|o  if  «  ;tO  • 

The  support  of  S\p{x)  =  p  *  (fiix)  is  larger  than  the  support  of  (t)(;c)  and  thus  larger  than  2.  This 
implies  that  S\p(x)  is  a  function  that  oscillates  between  its  zero-values  at  integer  points.  At  the 
finner  scale  2',  the  wavelet  transform  Wjjix)  therefore  has  local  maxima  at  non-intcgcr 
points.  Hence,  the  discrete  Dirac  is  an  example  of  discrete  signal  which  does  not  admit  an  inter- 
polation 5ip(x)  such  that  for  all  scales  2-'  (J  >  0),  W2,pix)  has  maxima  only  at  integer  points. 
Let  us  now  show  that  the  error  introduced  by  the  maxima  detection  can  be  viewed  as  a  high  fre- 
quency component  which  is  added  to  the  signal.  Let  h{x)  be  a  function  satisfying  the  condi- 
tions (87)  and  e{x)  =  Sip(x)- h(x).  We  can  choose  SiP(a:)  and  hix)  such  that  Uie  nonm  of 
the  error  ||e(;c)||  is  minimum.  The  function  e(x)  can  be  viewed  as  the  high  frequencies  of 
h(x)  that  cannot  be  written  as  a  convolution  of  some  function  with  (J)(a:).  Note  that  e{x)  is  zero 
at  every  integer  point.    Let   D  =    d„        _    be  a  discrete  signal  and  S\f{x)    be  the  function 

\.  J 

defined  by  equation  (89).  This  function  S\f(x)  is  an  interpolation  of  the  samples  of  D  and  we 
have 

Sifix)  =    2  ^1  5ip(x-n)  =    X  ^n  /^C^-")  +    Z   ^«  eU-«)  •  (91) 
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Since  h{x)  satisfies  the  conditions  (87),  we  proved  in  the  previous  paragraph  that  if  h(x)  could 

+00 

be  written  as  the  convolution  of  some  function  with  (5)(jt),  then  the  component     ^   dn  h{x-n) 

n  =—00 

would  have  a  finite-scale  wavelet  transform  whose  maxima  have  integer  abscissa.  The  second 
component     "^   d„   eix-n)  are  the  high  frequencies  of    ^   d„  h{x-n)  that  must  be  canceled 

in  order  to  obtain  a  function  that  can  be  written  /  *  ^{x).  These  high  frequencies  correspond  to 
the  error  introduced  when  detecting  the  maxima  at  integer  points. 

The  two-dimensional  extension  of  this  analysis  is  straight-forward  since  the  function 
<^{x,y)  can  be  decomposed  into:  ^{x,y)  =  <Sf{x)^(y).  Let  us  suppose  that  there  exists  a  function 
p(x)     which    satisfies    the    conditions    (87)    and    let    us    denote     R{x,y)  =  p{x)<;i{y).     Lei 


D  = 


dn. 


.  _„  and  let  us  define  the  function  f{x,y)&  L  (R^)  by 


f{x,y)=    X      I    d„^  9<.x-n)p{x-m)  .  (92) 

This  equation  yields 

Sxf{x,y)  =    Z      I.    d„^  SxR{x-n,y-m)  .  (93) 


From  the  conditions  (87),  we  derive  that  Sif{njn)  =  dn^.  Similarly  to  equation  (90),  we  can 
prove  that  Wj,/(x)  is  given  by: 


+00     400 


y^lf{x,y)=    X      I    Wlnnjn)  SxR{x-n,y-m)  .  (94) 

This  equation  implies  that  the  local  maxima  of  V/if{x  ,m )  along  x  have  integer  abscissa  and  are 
therefore  identical  to  the  discrete  maxima  of  Wl^f .  We  can  also  obtain  the  maxima  position  of 
^ifix,y)  along  x,  for  any  y  constant  from  the  maxima  for  y  integer.  The  same  observation 
can  be  made  along  the  vertical  direction  for  Wlf{n,y).  This  proves  the  extension  of  lemma  2 
that  was  stated  in  section  6.2. 


Appendix  5 
Projection  Operator  on  V 

In  this  appendix,  we  describe  more  precisely  the  projection  operators  Pr  defined  in  section 
3.2.  We  then  give  a  simple  implementation  of  this  operator  when  only  keeping  the  maxima  of  the 
wavelet  transform.  The  last  part  of  the  appendix  extends  these  results  in  two  dimensions. 
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Let  us  suppose  first  that  we  record  all  the  extrema  of  the  wavelet  transform 


W.Jix) 


;eZ- 


The  corresponding  set  T  regroups  all  sequences  of  functions 


gj{x)  such  that  ^/x)  is  in 


jeZ 


H  (R)  and  has  the  same  extrema  than  W2,f{x)  for  aU  y  e  Z.  These  conditions  impose  the 
value  of  gj{x)  on  a  discrete  set  of  points  and  since  gj{x)  must  not  have  any  other  extrema, 

-^ must  not  change  sign  between  two  extrema  of  W2jf  {x).  One  can  verify  that  these  pro- 
perties define  a  convex  set.  This  convex  set  is  however  not  closed  in  1  (L').  In  order  to  close 
r,  we  redefine  it  within  a  Hilbert  using  a  Sobolev  norm.  Let  us  recall  that  the  norm  of  a  function 
g{x)  in  the  Sobolev  space  H'(R)  is  defined  by: 


Let    1  (H  )    be  the   Hilbert  space  of  all  sequences  of  functions 
gj{x)s  HVR)  and 


[^My\, 


;eZ 


such   that 


The  value 


gjM 


J6Z 


defines  a  norm  in  this  Hilbert  space.  We  restrict  the  convex   T 


to  the  sequence  of  functions  in  1  (H  ).  In  this  Hilbert  space,  one  can  prove  that  the  convex  T  is 
closed. 


Let  us  now  define  the  operator  Pp  that  transforms  any  sequence     gj{x) 


;eZ 


l'(H>) 


into   the   closest  sequence 


hj{x) 


J;ez 


€  r,   with   respect   to   the   norm   of    I"(H  ).     Let 


tj {x)  =  hj{x)  -  gj {x ).  One  can  easily  prove  that  each  function  hj (x )  is  chosen  such  that 
lle>(^)llH-  =  Tl^>(^)l^^  +7l-^^|2^     isminimum. 


(95) 


Let  (xo.ao)  and  {x\,a\)  be  respectively  the  abscissa  and  amplitude  of  two  consecutive  extrema 
of  W^ifix).  Let  us  suppose  for  example  that  ao  <  ^i.  since  hj{x)  has  the  same  extrema  that 
Wj^/U),  wehave 


hj{XQ)   =   gjiXQ)   +  Zj{Xq)   =  Oq 

hjixi)  =  gjixO  +  tjixO  =  fli 

dhi(x)        dgiix)        deAx)   ^  ^   c  r  i 

— ^^-^  =     ^i,        +  — y-^  >  0    for  X  e  [xo  ,  a:  i] 


(96) 


dx 


dx 


dx 


dt.ix)  .-, 


The  problem  is  therefore  equivalent  to  the  minimization  of  J  \tj{x)\'^  dx  +   \  \ — ^ — \- dx 
with  the  three  constraints 
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Zj{xq)  =  aQ  -  gjixo) 
EjUi)  =  a  I  -  gj(xi) 
dt,(x)    ^   -dg,(x)     f  r  , 

^        -   — W~^  ^  ^  [Xo,Xi] 


(97) 


This  minimization  of  a  Sobolev  norm  with  convex  constraints  admits  a  unique  solution  which 
defines  £jix)  for  x  e  [j:o,j:i].  The  same  argument  can  be  used  between  -o°  and  the  first 
extrema  and  the  last  extrema  and  +«>  by  observing  that 

}}mjj{x)  =    \}m£jix)  =    0  .  (98) 

We  have  therefore  uniquely  defined  the  function  £j{x)  which  characterizes  the  projector  Pp. 
Let  us  suppose  that  f  {x)  b  H'(R).  Similarly  to  the  energy  conservation  equation  (7),  we 


can  show  that 


l/(;c)||^,  =  'Z\\w^f{x)\\^, 


(99) 


This  proves  that  if  we  restrict  the  space  V  only  to  the  wavelet  transfomis  of  the  functions 
f{x)€  H'(R),  then  it  is  included  in  l^(H').  One  can  also  show  that  the  reproducing  tccmel 
equations  (11)  define  an  orthogonal  projection  Pv  from  1  (H  )  on  this  restriction  of  V. 
Hence,  the  operator  P  =  Pr  o  Pv  restricted  to  l^(H')  is  the  composition  of  a  nonexpansivc  pro- 
jection on  a  closed  convex  set  and  an  orthogonal  projection.  This  guaranties  that  the  iteration  pro- 
cedure converges  weakly  to  some  element  at  the  intersection  of  F  and  V.  Note  that  we  do  not 
give  any  convergence  rate  neither  do  we  prove  that  the  intersection  is  reduced  to  a  unique  point. 

If  instead  of  recording  all  the  local  extrema  of  W2,f{x)  we  only  detect  the  points  where 
I  W2,f  {x)\  is  maximum,  the  corresponding  set  T  is  not  convex  any  more.  Indeed,  suppose  that 
ixo,ao)  and  {x\,ai)  are  two  consecutive  maxima  points  of  Wjjix)  with 


WjJixo)  =  ao  >  0   and    ^^/(j:,)  =  ai  >  0 


Let 


g/ix) 


jeZ 


and 


gjHx) 


JjeZ 


be  two  sequence  of  functions  in  F.  By  definition,  (xo.<3c 


and  {xi,ai)  are  maxima  of  gUx)  and  gHx).  Hence,  in  the  interval  [a:o  ,j:i],  both  g}{x) 
and  gHx)  have  a  positive  minima.  One  can  easily  build  two  such  functions  such  that  if 
gj\x)  =  l-(g/(x)  +  gjHx))  \hengUx)  has  also  a  maximum  in   ]a:o  ,  j:i[.  This  implies  that  the 


sequence 


gjHx) 


jeZ 


is  not  in  T  and  therefore  that  V  is  not  convex.  Since  F  is  not  convex, 


we  cannot  build  an  nonexpansive  projection.   However,  we  can  stiU  define  an  operator  Pp  thai 


transforms  any  sequence 


gjix) 


J€Z 


e  l^(H')   into  the  closest  sequence      hj{x) 


;6/- 


F 


Let  (xo,ao)  and  (xi,aO  be  respectively  the  abscissa  and  amplitude  of  two  con.sccutive  maxim. 
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points  of  \W2,fix)\.  The  function  Zj(x)  =  hj(x)-  gj{x)  is  defined  by  the  minimization  of 

]  |e,(;c)|2^  +I\^l^i2^  ,  (100) 

with  the  constraints 

lEjixo)  =  ao  -  gjixo) 

\zjix,)  =  ax-  gjixO  (101) 

\Ej{x)  + gj(x)\  =  \hjix)\    must  have  no  maxima  for  xe  ]xo,j:i[  . 

One  can  show  that  this  minimization  problem  admits  a  solution  although  it  might  not  be  unique. 

Let  us  now  describe  the  discrete  implementation  of  the  operator  Pr  for  the  wavelet 
transform  maxima.  In  this  implementation,  we  simplify  slightly  the  definition  of  the  projector 
Pr  in  order  to  speed-up  the  computations.  We  compute  the  function  Ej  (x  )  which  minimizes 
(100)  given  the  constraints  of  equations  (101)  and  then  verify  that  \hjix)\  has  no  maxima  for 
x  e  ];co  ,x\[.  If  \hj{x)\  has  some  some  spurious  maxima,  we  remove  these  with  a  simple  cut- 
ing  procedure.  In  that  case,  the  function  £j(x)  is  not  optimal.  Within  our  discrete  model,  we 
showed  in  appendix  4  that  W2jf  (x)  can  be  written 

W:!,f{x)  =    X   ^2-/(«)  Sipix-n)  .  (102) 

f 

We  can  therefore  restrict  the  set  V  to  the  sequence  of  functions     gj  (x) 
j  €  Z,  gj  {x )  can  be  written 

gjix)  =     I   ^;(«)  Sipix-n)  .  (103) 


.,  such  that  for  any 


^  which  has  the  same 


This  is  cleariy  equivalent  to  specifying  the  discrete  sequences    gjin) 
discrete  maxima  than  W^f .  The  function  e,  (x)  can  also  be  written 

^j(x)  =    £  e;(n)Sip(x-/j)    . 

Let  us  add  this  constraint  to  the  minimization  problem  on  tjix).  Let  (/io,ao)  and  (Aii,ai)  be 
two  consecutive  maxima  of  W  2,  fix).  The  function  tj  (x )  should  minimize 

J  I  e,(;c)  |2  ^  +  ]  I  ^^  1 2  dx    .  (104) 


with  the  constraints 


tjino)  =  ao-  gjino) 

ejinO  =  ai  -  gjini)  ■  ^^^^^ 
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Let  us  suppose  that  the  function  Sipix)  is  given  by  equation  (37).  The  Euler  equations  associate 
to  this  minimization  problem  yield: 


Ejino)  =  ao  -  gj(no) 
^iej(n+l)-2£j{n)  + 
Ejini)  =  a\  -  gjini) 


^  {ej{n+l)-2£jin)  +  ej(n-\))-ejin)  =  0     for   no<n<ni 


We  compute  Zj(n)  for  no<  n  <  ni  by  solving  this  tridiagonal  linear  system  which  can  be  done 

with  standard  methods  in  0(na-n\).  If  the  resulting  discrete  sequence  \hj(n)    „^<„_  has  some 

spurious  discrete  maxima,  we  remove  these  by  setting  the  signal  to  a  constant  on  the  righi 
domain.  We  find  the  appropriate  rib  and  rie  such  that  the  spurious  maximum  belongs  to  [nt,  ,  n^  ] 

and  set  all  samples  g{n)  for  n  e  [rib  ,n«]  to         ^   ~        ^    .  This  maxima  removal  procedure 

has  also  a  complexity  (9(/io-«i).  If  the  original  discrete  signal  D  has  A'  non-zero  samples, 
each  discrete  signal  W^/  has  also  A'  non-zero  samples.  The  computation  of  Zjix)  thus 
requires  0{N)  computations.  Since  there  are  at  most  log(,N)  different  scales,  the  implementa- 
tion of  the  projector  Pf*  has  a  complexity  0{N  [og(N)). 

In  two-dimensions,  the  implementation  of  the  non-linear  projection  operator  Pj^  is  based 
on  the  one-dimensional  projector  that  we  just  described.  The  maxima  of  the  wavelet  transform 
are  detected  along  the  rows  of  Wj/^f  and  the  columns  of  W^-''/ .  At  a  given  scale  2-' ,  for  the 
horizontal  orientation,  the  projector  Pf^  associates  to  a  discrete  two-dimensional  signal  ¥/ 
another  discrete  signal  Z/  which  has  the  same  maxima  than  Wl-^f  along  each  of  its  rows.  The 
continuous  model  behind  this  discrete  operator  is  a  simple  two-dimensional  extension  of  the 
one-dimensional  model.  Each  row  of  Yj^  can  be  viewed  as  a  one-dimensional  discrete  signal 
whose  maxima  must  be  adapted  in  order  to  match  the  maxima  of  the  corresponding  row  of 
Wl'^f .  The  problem  is  therefore  reduced  to  a  one-dimensional  processing  and  we  can  directly 
use  the  discrete  operator  Pf^  defined  for  one-dimensional  signals.  In  two  dimensions,  the  imple- 
mentation of  the  operator  Pf*  thus  consists  of  applying  the  one-dimensional  algorithm  along  each 
row  of  Wl'"^/  and  columns  of  W^-^f  for  all  scales  2-'.  One  can  easily  check  that  the  complex- 
ity is  still  O  (N  log  {N )  where  N  is  the  total  number  of  pixels  of  the  original  image  D. 


Appendix  6 

A  Particular  Qass  of  Two-Dimensional  Dyadic  Wavelets 

In  this  appendix  we  characterize  the  two-dimensional  wavelets  used  for  numerical  computa- 
tions. The  two  wavelets  H'\x,y)  =  \\f{x)%(y)  and  H'^U ,>- )  =  ^(x )  vCy )  are  defined  from  a  one- 
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dimensional  wavelet  \|/(x)  satisfying  the  conditions  of  appendix  2.  The  Fourier  transform  yCco) 
thus  satisfies 

\iIr(2co)  =  G((0)<^((0)    with   (^(0))  =  n//(2-''(0)  .and  (107) 

|G((0)|2  +   |//(co)|2  =  1  .  (108) 

One  can  easily  show  that  the  Fourier  transforms  4^'(C0:t,C0y)  and  4'^(C0i,C0j,)  satisfy  the  condi- 
tion (61)  if  and  only  there  exists  a  function  0(j:  ,y )  whose  Fourier  transfonm  satisfies: 

\4>((o,,0iy)\^  -   |6(2cO;,,2cOy)|2  =  |^'(2Q);,,2cOy)l'  +   |T'(2(0;,,2cOy)|2  ,  (109) 

and  lim         |<t>(C0;c,(av)|  =  1  . 

From  equations  (107)  and  (108),  one  can  prove  that  (109)  holds  for 

if  and  only  if: 

|G(coJ|2  \mx)\^  lC(2a),)|2  +  |^(2co,)|2  |G(co,)|2  |(^(a),)|2  = 

(1  -  |//((0,)|2  |//((0y)|2)  |(^((0J|2  |$(C0,)|2   .  (110) 

The  Fourier  transform  ^(co)  must  therefore  satisfy 

C(2co)  =  L((o)$(co)    with    |L(a))|2  =    ^+1^(^)1'    .  (Ill) 

The  function  L(co)  can  be  interpreted  as  the  transfer  function  of  a  discrete  filter  L.  Since  (^{x) 
is  a  smoothing  function,  one  can  derive  that  ^(co)  is  the  Fourier  transform  of  a  smoothing  func- 
tion t,ix).  The  Fourier  transform  of  the  functions  <t>(x,y)  ,  ^\x,y)  ,  and  ^^(x,y)  are  thus 
given  by: 

O(c0;,,t0y)  =  n//(2-''(Ox)//(2-''cOy)   , 
4''(2(0x,2cOy)  =  G((0;t)L(a)y)O(m^,0)y)    and 

4'2(2C0i,2C0y)  =  L(Mx)G((05,)O((0;c.0>y)    . 


^  whose  transfer  function  is  the  func- 

neZ 


Table  4  gives  the  impulse  response  of  a  filter  L  = 

tion  L((o)  computed  with  //(co)  given  by  equation  (83).  This  filter  corresponds  to  wavelets 
^^{x,y)  and  ^'^{x,y)  built  from  the  one-dimensional  cubic  spline  wavelet.  The  filter  L  is  sym- 
metrical aroimd  zero. 
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n 

In 

0 

0.79113 

1 

0.06623 

2 

0.03118 

3 

0.00727 

4 

0.00003 

Table  4:  First  five  coefficients  of  the  impulse  response  of  a  filter  L    corresponding  to  the  cubic 
spline  wavelet. 

Table  5  gives  the  impulse  response  of  a  filter  L  corresponding  to  the  function  //((o) 
given  by  equation  (84).  This  filter  corresponds  to  wavelets  ^^{x,y)  and  4''(j:,.v)  built  from  the 
one-dimensional  Haar  waveleL  It  is  also  is  symmetrical  around  zero. 


n 

In 

0 

1 
2 
3 

0.8598 

0.0729 

-0.0031 

0.0002 

Table  5:  First  four  coefficients  of  the  impulse  response  of  a  filter  L   corresponding  to  the  Haar 
wavelet. 


Appendix  7 
Fast  Wavelet  Algorithms  for  Two-Dimensional  Signals 

We  describe  an  algorithm  for  computing  a  discrete  wavelet  transform  and  the  inverse 
wavelet  transform  algorithm  for  images.  We  suppose  that  the  two  wavelets  4^'U.y)  and 
T^Cx,}')  are  characterized  by  the  three  discrete  filters  H  ,  G  and  L  mentioned  in  appendix  6.  Wc 
denote  by  Hp  ,  Gp  and  Lp  the  discrete  filters  obtained  by  putting  1p  -  1  zeros  between  each 
coefficients  of  respectively  the  filters  H  ,  G  and  L  .  The  transfer  function  of  these  filters  is 
H(2Poi),  GOP  (a)  and  L(2P(d).  We  also  denote  by  Hp  ,  Gp  and  Lp  the  filters  whose  transfer 


function  are  respectively  H{2P(m)  ,   G(2P(o)   and  L(2Poi)  .    We  denote  by   A  *  (H ,L)    the 
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separable  convolution  respectively  of  the  rows  and  columns  of  the  image   A    with  the  one- 
dimensional  filters  H  and  L. 

The  following  algorithm  computes  the  two-dimensional  discrete  wavelet  transform  of  an 
image  Sif .  At  each  scale  2-' ,  the  algorithm  decomposes  5^.,/  into  S^f  .'W]j'^f  and  V^]j'^f  . 

;  =  0 

w/zi/e  (}  >  J) 

Wi-^J  =  Sif  *(Gjd.j) 

Wl-if  =Sif  HLj,Gj) 

Si4  =  Sif  *{HjMj) 

j=j  +  } 
end  of  while 

The  proof  of  this  algorithm  is  based  on  the  properties  of  the  wavelets  H'^(x  ,y)  and  *f-(j:,>') 
described  in  appendix  6.  If  the  original  image    Sif{n,m)  has  A'   non-zero  pixels,  the 

complexity  of  the  algorithm  is  O  {N  log  (N )).  In  order  to  avoid  border  problems,  one  can  extend 
the  signal  by  symmetry  around  the  first  and  last  rows  and  the  first  and  last  column. 

As  in  the  one-dimensional  case,  the  reconstruction  algorithm  computes  Sif  by  recon- 
structing at  each  scale  2^  the  signal  S^.J  from  S^f  ,  Wj//  and  Wl-^f .  The  complexity  of  this 
reconstruction  algorithm  is  also  O  (N  log  (N )). 

J  =  J 

while  (j  >  0) 

Si-J  =  V^Vf  *  (G;-iA-i)  +  W2^/  *  (L,_,.G,-,)  +  Sif  *  (Hj-iJij-0 

end  of  while 
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